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SUMMARY 

^ * 

This  report  describes  the  research -ctm4u£Lted  under  the 

Office  of  Naval  Research  Contract  No.  ^00*14-77 -C03T^to  in- 
vestigate algorithms  for  locating  simultaneously  multiple 
signal  sources/ jammers  using  linear  antenna  arrays  and  digi- 
tal signal  processing.  Specifically,  the  signals  impinging 
on  each  of  the  antenna  elements  are  converted  to  base  band 
and  digitized.  The  digital  samples  from  all  the  antenna  ele- 
ments are  fed  to  a computer  which  carries  out  direction  find- 
ing, and  later  beam  forming  and  signal  extracti on , (Fi gure  S-l) 

^ The  problem  of  multiple  signal  source  location  with  a 
linear  array  -ca-n  -ke  formul  ated.  ars^ol  1 ows  . Consider  the 
linear  array  with  n isotropic  antenn\  el ements  spaced  at  in- 
tervals of  d units.  The  array  is  operated  in  a narrowband 
mode  at  a center  frequency  whose  wavelength  is  y.  There  are 
m signal  sources/ jammers  0,,  Jn,  ...»  Jj.  located  with  inclina- 


tions of  0. 


0, 


0. 


1*  u2’  - ’ -m 
and  complex  envelope  amplitudes  s 


1*  2 nft 

with  the  antennia  normal  (Figure  S-2) 


1*  2’ 


m 


Let  the 


complex  envelope  amplitude  of  the  antenn^  internal  noise 

. Assume  E { g ^ } 


generated  by  the  ith  antenna  element  be 
<$ i j g 2 and  Etgts.}  = 
simplicity  E { s . } = 0 and  A = ||E{s* 


0,  E{g*9j} 


Iso  assume  for 


Let  z..  be  the 


combined  signal  output  of  antenna  i due  to  the  signal  sources 

and  internal  antenna  noise  and  C = ||E{z*z.}||.  Given  C the 

i J 

problem  is  to  evaluate  the  number  of  signal  sources  and  their 
angular  locations.  (^nce  this  problem  is  solved,  it  is  a 
simple  matter  to  extend  to  the  case  where  two  mutually  per- 
pendicular linear  arrays  are  used  to  determine  azimuth  and 
elevation  angles  of  signal  sources^. 


t\1 


KM 


Figure  S-l.  Block  diagram  of  radar  receiver  system  with  digital  signal  processing 
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In  order  to  solve  the  problem  of  multiple  source  loca- 
tions, a theoretical  basis  is  developed  using  the  properties 
of  Vandermonde  determinants  and  vectors  and  several  algorithms 
are  formulated  utilizing  this  basis.  Basically  these  algorithms 
consist  of  three  computational  phases.  In  the  first  phase, 
eigenvalues/eigenvectors  of  a complex  covariance  matrix 
formed  from  the  sampled  output  signals  of  the  elements  of 
the  linear  array  are  computed.  In  the  second  phase,  a poly- 
nomial is  formed  from  the  eigenvectors  computed  in  the  first 
phase  and  in  the  final  phase  the  polynomial  is  solved  on  the 
unit  circle  to  find  the  angular  locations  of  the  signal  sources. 

The  proposed  algorithms  are  applicable  to  cases  where  the 
source  signal  covariance  matrix  A is  singular  as  well  as  non- 
singular. The  source  signal  covariance  matrix  is  singular, 
for  instance,  in  the  presence  of  the  multipath  phenomenon. 

Using  a first  order  perturbation  analysis,  the  effects  of 
errors  due  to  quantization  and  unidentical  variances  and  non- 
zero covariances  of  antenna  internal  noise  signals  are  dis- 
cussed and  a method  is  suggested  to  minimize  these  effects. 

Also  an  upper  bound  is  derived  on  the  computational  errors 
in  the  estimates  of  angular  locations. 

A computational  analysis  is  presented  to  show  that  the 

•5 

proposed  algorithms  have  a computational  complexity  of  0(n  ) 
multiplication.  The  most  complex  and  generalized  version 
of  the  proposed  algorithms  requires  approximately  175  million 
multiplications  for  n = 100  and  around  50  million  multiplica- 
tions for  n = 60.  Thus  if  angular  location  of  targets  has 
to  be  performed  every  10  milliseconds,  we  require  a proces- 
sing rate  of  approximately  17500  and  5000  million  multiplica- 
tions/second for  n = 100  and  60  respectively.  (These  rates 
have  to  be  multiplied  by  2 if  elevation  and  azimuth  angles 
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of  the  sources  are  required).  In  view  of  the  present  tech- 
nological trend  in  computing,  it  is  feasible  to  implement 
the  algorithms  for  real  time  angular  location  on  high  speed 
parallel  processing  computers  within  the  next  decade  when 
n is  in  the  range  10  - 100. 

Simulations  are  performed  to  investigate  the  behavior 
of  the  algorithms  in  terms  of  their  accuracy,  resolution  and 
convergence.  The  signal  sources/ jammers  and  the  internal 
noises  are  simulated  by  generating  independent  normally 
distributed  random  numbers.  Specifically  experiments  were 
conducted  to  illustrate  that  in  general  as  the  number  of 
samples  and/or  the  number  of  antenna  elements  are  increased, 
accuracy  of  angular  estimates  improves.  Also,  the  experi- 
ments illustrate  that  higher  number  of  elements  and/or 
higher  signal  power  levels  result  in  better  angular  reso- 
lution of  signal  sources.  Preliminary  experimental  results 
indicate  that  for  four  antenna  elements,  for  widely  separated 
signal  sources,  an  accuracy  of  approximately  10  mllliradians 
can  be  expected  after  30-40  samples.  Also  for  four  elements, 
a resolution  of  approximately  7 degrees  may  be  expected  with 
the  source  signal  to  antenna  internal  noise  ratios  around  10 


In  order  to  assess  the  effectiveness  of  the  algorithms 
under  real  life  conditions,  actual  test  data  tapes  were 
obtained  from  the  Naval  Research  Laboratory  and  the  algorithms 
were  applied  to  process  data.  These  tapes  contained  covariance 
matrices  obtained  from  a 4-element  linear  array.  The  number 
of  sources  varied  from  0 to  10  and  an  anal og-to-di gi tal  con- 
version scheme  with  512  samples  was  used.  To  compute  the 
covariance  matrices,  1024  samples  were  used. 


The  covariance  matrices  for  the  cases  where  the  number 


of  sources  Is  from  1 to  3 were  processed.  The  results  are 
encouraging  for  one  and  two  sources  (except  for  certain  two 
source  cases)  and  the  angular  errors  are  below  5 degrees. 
However  for  three  sources,  the  algorithm  performs  less  satis- 
factorily in  terms  of  resolution  coalescing  two  or  three 
sources  which  are  separated  approximately  by  10  degrees  into 
one.  This  is  mostly  due  to  the  reasons  that  the  data  is  ac- 
curate only  up  to  the  second  decimal  digit  and  that  there 
are  only  four  antenna  elements. 

The  proposed  entirely  digital  approach  to  multiple  source 
location  (as  opposed  to  the  hybrid  analog-digital  approach  de- 
veloped by  Berni  and  Swarner  of  the  Ohio  State  University*) 
should  be  of  interest  to  radar  system  engineers  who  are 
faced  with  the  problem  of  detecting  targets  in  the  presence 
of  heavy  jamming  and  to  electronic  navigation  specialists. 


Swarner,  W.G.  and  Berni,  A.J.,  An  Adaptive  Antenna  Array 
for  Angle  of  Arrival  Estimation  System  for  Sensor  Com- 
munications, Report  No.  ESL3435-2,  Electro  Science  Laboratory, 
Department  of  Electrical  Engineering,  The  Ohio  State  University, 
Columbus,  Ohio.  Also  see  Berni,  A.O.,  "Angle  of  Arrival  Esti- 
mation Using  An  Adaptive  Antenna  Array",  IEEE  Trans.  Aero- 
space and  Electronic  Systems,  Vol . AES-11,  No.  2,  March, 
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1.  INTRODUCTION 


This  report  describes  the  results  of  a research  study 
conducted  to  find  and  explore  algorithms  for  angular  loca- 
tion of  multiple  signal  sources  employing  digital  signal 
processing  and  linear  antenna  arrays.  A theoretical  basis 
is  presented  using  the  properties  of  Vandermonde  determinants 
and  vectors  and  several  algorithms  are  developed  using  the 
basis  for  simultaneous  angular  location  of  multiple  signal 
sources.  For  the  algorithms  it  is  assumed  that  the  internal 
antenna  noise  signals  in  the  individual  elements  of  the  linear 
array  have  identical  variances  and  are  uncorrel ated . Basically 
the  algorithms  consist  of  determining  the  eigenvalues/eigen- 
vectors of  a complex  covariance  matrix  obtained  from  the 
sampled  output  signals  of  the  linear  antenna  array  and  solv- 
ing a polynomial  equation  whose  coefficients  are  obtained 
from  the  eigenvectors. 

The  proposed  algorithms  are  applicable  to  cases  where 
the  source  signal  covariance  matrix  is  singular  as  well  as 
nonsingular.  (The  source  signal  covariance  matrix  is  singu- 
lar when  the  multipath  phenomenon  is  present).  The  effects 
of  errors  due  to  quantization  and  unidentical  variances  and 
nonzero  covariances  of  antenna  internal  noise  signals  are 
discussed  using  a first  order  perturbation  analysis.  A 
computational  analysis  is  presented  to  show  that  the  proposed 
algorithms  have  a complexity  of  0(n  ) multiplications  where 
n is  the  number  of  elements  in  the  linear  array  and  that  it 
is  feasible  to  implement  them  on  high  speed  parallel  proces- 
sing computers  for  real  time  angular  location  In  the  next 
decade.  Simulations  are  performed  to  investigate  the  be- 
havior of  the  proposed  algorithms  in  terms  of  their  accuracy, 
convergence  and  angular  resolution.  Actual  test  data  was 
obtained  from  the  Naval  Research  Laboratory  to  test  the  al- 
gorithms and  the  results  are  reported.  The  report  is  concluded 
with  suggestions  for  future  research  in  this  area. 
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2.  GENERAL  THEORY  OF  MULTIPLE  SOURCE  LOCATION 


? 

The  problem  of  simultaneous  angular  location  of  multiple 
sources  with  a linear  antenna  array  can  be  formulated  as  fol- 
lows. Consider  the  linear  array  with  n isotropic  antenna  ele- 
ments spaced  at  intervals  of  d units.  The  array  is  operated 
in  a narrow  band  mode  at  a center  frequency  whose  wavelength 
is  v.  There  are  m signal  sources/ jammers  J^,  J2>  ....  Jm 
located  with  inclinations  of  0^,  02»  ...»  0m  with  the  antenna 
normal  (Fig.  2-1)  and  complex  envelope  amplitudes  s^  s2>  ...» 
sm.  Let  the  complex  envelope  amplitude  of  the  internal  noise 
generated  by  the  ith  antenna  element  be  g^  The  following 
first  and  second  order  statistics  are  assumed  for  the  signal 
and  antenna  noise  waveforms: 

E (g.)  = 0 

E l9)9j>  - S,^2 
<i  » 

e {g^}  = o 

Assume  for  simplicity  E { s ^ } = 0 and  let  A be  the  source  signal 
covariance  matrix: 

A - I la-jl  | = I I E( s*s j > | | 

Associate  a direction  vector  with  J..  as  follows: 


I 


X 

2 


Lfs(=&: 

*1 = 


exp  { j ( 2 tt  d / v ) sin  9-} 
exp  {j  (4TTd/v)  sin  0 ^ > 


exp  (j(n-l)  ( 2 -rr  d / v ) sin  0 . } 


Let  zT  = (Zp  z2>  ...»  zm)  where  z..  is  the  combined  signal 
output  of  antenna  i due  to  the  signal  sources  and  internal 

antenna  noise  and  gT  = (g^,  g^ gm).  Then  z can  be 

written  as: 


■ 9 + l S1 
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and  the  covariance  matrix  of  the  combined  signal  vector 
as: 

C = ||EU*Zj>|| 

■ g2!  + I U) 

i . j 

Given  C the  problem  is  to  evaluate  the  number  of  signal 
sources  and  their  angular  locations  using  the  relation 
(1). 

To  solve  this  problem  we  make  use  of  certain  remark- 

2 

able  properties  of  the  vectors  of  the  form  (1,  a,  a , ...» 
k- 1 

a ) referred  to  as  Vandermonde  vectors  in  the  following. 
(Note  the  direction  vectors  are  Vandermonde  vectors.)  We 
denote  such  vectors  by  V^(a)  or  simply  V(a).  We  need  the 
following  properties  of  Vandermonde  vectors  for  our  de- 
velopment. 

Lemma  1:  Vectors  V^a.),  i = 1,  2,  ....  I,  are  linearly 
independent  for  k ^ l and  a.  f a.,  i f j. 

The  lemma  is  a consequence  of  the  properties  of  Vander- 
monde determinants  (NERI  67).  This  simple  lemma  has  certain 
interesting  consequences.  Let  {Xj,  x2,  ...»  x^}  be  a set 
of  linearly  independent  vectors  with  k components,  k =1 
and  S be  the  space  by  these  vectors.  Assume  that  S possesses 
a basis  V which  consists  entirely  of  Vandermonde  vectors, 

V = V^(a2),  Such  a basis  will  be  re- 

ferred to  as  Vandermonde  basis  and  it  can  be  seen  from  Lemma 
1 that  if  a Vandermonde  basis  exists,  it  is  unique. 
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Theorem  1:  Let  (Xj»  Xg Xj^  be  a set  of  linearly  in- 

dependent k-tuples  with  Vandermonde  basis  {V^c^),  Vk(a2), 
. . . , V k ( ot ^ ) } , k > i.  Let 

F " (X  ^ |x  2 I * * * lx  4! 


f i 

|i 

li 


Construct  polynomials: 

G ( x ) * wj  (W+W)“1WtVk_1(x)  - xk-1 
gt(x)=  w!  (w+w)“1w+vk_1(x)  - x1"1, 

1 = 1,2 k-i 

Then  f (x)  - g.c.d.{G(x),  g1(x) gk-1(x)}  possess  aj, 

«2 a£  as  its  only  roots. 

Proof:  The  existence  of  Vandermonde  basis  Implies  the 
existence  of  i ( * + 1 ) -tupl es  (a,  y^  y? yfc)  satisfy- 

ing the  following  k nonlinear  equations: 
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T 


vk(a)  - l y,x, 

i 

i.e.  V k (a ) = Fy 


whore  y ^ (y x , y2. 


y^)1^.  Rewrite  (2)  as: 


(2) 


Wy  =>  vk_1(o) 

wjy  = ak~l  (3) 

Since  a consistent  solution  y exists,  y can  be  written  as  y = 
(W!W)  1Wl’vk  l(a)  and  hence  (3)  becomes 

G (a)  = wJ(W+W)"1W+Vk_1(a)  - ctk-1 

G(a)  possesses  as  its  roots  and  ( k - l - 1 ) 

extraneous  roots.  To  eliminate  the  extraneous  roots,  we  note 
for  the  desired  roots: 

Wy  = V k_  j (a ) 

i.e.,  W(W  + W)"1W+Vk_1(a)  = V|(_1(a) 

(W(WtW)"1W+  - I)  Vkl(a)  = 0 

from  which  one  can  obtain  polynomials  g^(a),  g2(a),  .... 
9k_i(a)-  Thus  f(x)  - gcd{G(x) , g^x)}  has  as  its  roots 
ctj,  a^,  ....  a^.  It  can  be  seen  that  if  a*  is  a root  of 
f(x),  then  y can  be  found  so  that  Fy  = Vk(ct*)  and  in  view 
of  Lemma  1,  a*c{ot^,  •••»  Thus  the  only  roots  of 

f(x)  are  oi^,  , ...»  a^. 
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We  refer  to  G(x)  as  the  principal  polynomial  and  g..(x) 
as  auxiliary  polynomials  in  the  sequel. 

Corollary  1:  For  k = l + 1,  (x)  = 0 and  hence  f(x)  = G(x). 

Hence  a set  of  l.i.  ( £+1 ) -tupl es  {X^  X,,,  X^}  always  pos- 

sess a unique  Vandermonde  basis  (for  f(x)  is  of  degree  l) . 

The  corollary  is  interesting  since  it  uniquely  character- 
izes Si  vectors  by  means  of  an  A-tuple  (cij,  a 2 a^). 

Corollary  2:  If  vectors  X^,  X2>  X^  are  orthonormal, 

then 

G ( x ) = wJw+Vk_1(x)  - (1-X)xk_1 

g^x)-  w{(I  + (l-X)'1w*w[]W+Vk_1(x)  - x1'1 
2 

where  \ = ( w k | . 

This  follows  by  noting  W+W  = (I-w*wT)  and  applying 
Sherman-Morri son  formula  (DAHL  74)  to  obtain  its  inverse. 

The  implication  of  Corollary  2 is  that  it  is  computational- 
ly simpler  to  obtain  f(x)  when  the  vectors  Xj,  X2>  ....  X^ 
are  orthonormal . 
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2.1  Multiple  Source  location  for  Non-Singular  Source 
Signal  Covariance  Matrix 


The  problem  of  multiple  source  location  when  the  source 

signal  covariance  matrix  A is  non -singular  is  simpler  and 

will  be  considered  first.  We  show  that  there  are  exactly 

2 

m eigenvalues  of  C which  are  greater  than  g and  that  for 
their  corresponding  m eigenvectors,  the  in  source  direction 
vectors  form  a basis.  Then  we  apply  Theorem  1 to  obtain 
the  direction  (Vandermonde)  vectors  and  consequently  angular 
1 ocations . 


Consider  the  combined  signal  covariance  matrix  C which 
can  be  written  as: 


• » J 

= g2I  + 5A5+ 


where  H = [ £ ^ | £2 I S3  I • • • I • To  show  that  C has  m eigen- 
values greater  than  g2 , we  show  that  the  received  signal 

t 

covariance  matrix  B = HAH  has  m positive  eigenvalues  and 
n-m  zero  eigenvalues.  We  also  show  that  the  eigenvectors 
of  B (and  hence  C)  can  be  expressed  as  Y z.£.  where  z^  = 

li 

i 

(Zj,  z2,  ....  zfl))  is  an  eigenvector  of  AH  S=  AA. 

Let  n = l zi^i*  *n  or£*er  that  n is  an  eigenvector 
i 

of  B,  we  should  have: 
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Bq  3 An 


*J* 

I.e.,  HAH  Hz  = AHz 
HAAz  = AHz 
H ( AA- A I ) z - 0 
<==«>  (AA-AI)z  = 0 

Thus  the  eigenvalues  of  B are  the  same  as  those  of  AA. 

It  may  be  noted  that  A is  Gramian  of  H (GANT  60)  and  has  a 
rank  m and  since  A also  has  a rank  m (for  it  is  non-singular) 
it  follows  that  B has  m nonzero  eigenvalues.  It  also  follows 
that  these  eigenvalues  are  positive  since  A and  A * are  both 
positive  definite  ( W I L K 66). 

Let  c .,  ....  r be  the  eigenvectors  of  B that  have 

eigenvalues  greater  than  g'.  Since  these  eigenvectors  can 
be  chosen  to  be  orthonormal  (for  B is  Hermitian),  we  apply 
Corollary  2 to  obtain  the  direction  (Vandermonde)  vectors. 
(Note  n > m for  the  theorem  and  corollary  to  be  applicable). 
We  first  obtain  G(x)  and  then  solve  the  polynomial  equation 
to  obtain  its  roots.  In  this  case  since  the  roots  happen  to 
be  of  the  form  exp  (-j(2ird/v)  sin  0),  we  restrict  our  search 
for  the  roots  to  points  on  the  unit  circle.  This  simplifies 
matters  considerably.  (It  is  our  experience  that  the  extran- 
eous roots  of  the  principal  polynomial  G(x)  rarely  lie  on  the 
unit  circle.  Thus  the  auxiliary  polynomials  and  their  roots 
are  not  usually  computed.  Later  on  we  consider  the  entire 
set  of  polynomials  (G(x),  g^(x)}  as  an  overdetermined  non- 
linear system  (in  one  variable)  and  compute  least  squared 
error  solutions.) 

Ml 


10 


Example  1:  Five  antenna  elements  are  assumed  spaced  at  half 
wavelengths . To  simplify  matters  the  covariance  matrix  of 
antenna  internal  noise  is  assumed  to  be  an  identity  matrix. 
Three  sources  are  simulated  at  inclinations  of  0,  .In  and 
-.In  to  the  antenna  normal.  The  source  signals  are  assumed 
to  be  samples  from  white  Gaussian  noise  and  generated  by 
using  Box-Muller's  transformation  (HAUL  74).  The  signal 
variances  are  assumed  to  be  5 power  units.  After  10  samples, 
the  received  signal  covariance  matrix  is  evaluated  and  an 
identity  matrix  is  added  to  obtain  the  combined  signal  co- 
variance  matrix.  The  eigenvalues  of  the  combined  signal 
covariance  matrix  are  32.8,  14.34,  8.67,  1.0,  1.0  and  the 
eigenvectors  corresponding  to  eigenvalues  having  a value 
greater  than  1 are: 


02368  +j. 37984 


-.09574  + j. 11358 


-.26674  -j. 33110 


-.34044  -j. 56700 


-.25265  -j. 38871 


.42832  -j. 38035' 


-.  53024  -j. 37964 
-.41012  -j. 17065 


-.17255  +j. 06464 


,02439  +j. 12137 


-.62973  + j. 03864 
.16043  - j. 01170 
.40649  - j .0096 
-.00581  + j. 04302 
-.71746  + j . 10033 


The  principal  polynomial  G(x)  is  computed  as: 

G ( x ) = .19458  - . 16938x  - .10721x2 
+ . 3 2 6 9 5 x 3 - .24493x4 

There  are  three  roots  e3  0 a 0,  .97081  , 5.31238,  for  which 
G(x)  vanishes  and  the  three  angular  locations  are  obtained 
with  seven  digit  accuracy. 

Now  we  consider  the  situation  where  there  is  'noise' 
present  in  the  system  of  non-linear  equations  due  to  such 
factors  as  receiver  noise  and  the  non-whiteness  of  antenna 
noise.  In  this  case  the  desired  roots  of  G(x)  and  g^(x)  do 
not  exactly  lie  on  the  unit  circle  and  we  seek  to  find  the 
least  squared  error  solution  satisfying  the  system  of  non- 
linear equations.  We  write  the  system  of  nonlinear  equa- 
tions as  the  following  matrix  equation: 

FVr(x)  = 0 
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We  minimize  <{>  = Vk(x)  F 1 FVk(x)  subject  to  the  constraint  that 
x lies  on  the  unit  circle.  Let  x = and  F^F  = i l^i  ke^°llcl  I • 

Then  <j>  = £ ^exptjlO  . ^-(i-k)OJ } . To  find  minima,  we 

i ,k 

set  34./30  = 0 and  32<j>/302>0  to  obtain: 

4>'(0)  = 2 l qik(i-k)sin(e.k-(i-k)0)  = 0 and 
i>k 

(j>"(0)  = -2  l q-k(i-k)2cos(0 -k-(i-k)0)>0. 
i >k 

In  next  section,  we  explain  a method  which  does  not 
entirely  rely  on  the  values  of  eigenvalues. 


2.2  Multiple  Source  Location  for  Singular  Source 
Signal  Covariance  Matrix 


The  problem  of  multiple  source  location  becomes  complex 
if  the  source  signal  covariance  matrix  A happens  to  be  singu- 
lar. The  singularity  of  A may  be  because  of  the  multipath 
phenomenon  or  the  deliberate  enemy  control  and  transmission 
of  jamming  or  signal  waveforms.  The  complexity  of  the  prob- 
lem is  due  to  the  fact  that  there  will  be  only  m'<m  eigen- 
values  of  C which  are  greater  than  g and  the  corresponding 
m'  eigenvectors  do  not  form  a basis  for  the  m source  direc- 
tion vectors;  here  m*  is  the  rank  of  A.  Our  approach  to  the 
problem  is  to  generate  eigenvectors  for  different  covariance 
matrices  obtained  by  considering  n,n+l,...  antenna  elements 
and  select  independent  vectors  from  these  sets  of  eigenvect- 
ors. These  vectors  form  a basis  for  the  source  direction 
Victors  and  hence  Theorem  1 can  be  applied  to  obtain  angular 
locations.  It  should  be  noted  that  this  approach  is  appli- 
cable when  A is  non-singular  and  the  eigenvalues  of  AA  are 

2 

small  and  comparable  to  g , thus  making  it  difficult  to  apply 
the  algorithms  suggested  in  the  preceding  section. 


Consider: 


C = g2I  + SAS+ 


We  show  that  C has  in'  eigenvalues  greater  than  g where  m'  = 
rank(A)  by  proving  that  B = HAH  has  m1  postive  eigenvalues 
and  n-m'  zero  eigenvalues.  As  before,  we  can  prove  that 
£ z.£.  is  an  eigenvector  of  B with  an  eigenvalue  of  X if 
1 T 

( z j , » ....  zm)  is  an  eigenvector  of  AA  with  an  eigenvalue 

of  X.  Since  A is  of  rank  m'  and  A has  a rank  m,  there  are 
m'  non-zero  positive  eigenvalues  of  AA  and  consequently  m' 
positive  eigenvalues  of  B.  Also,  one  may  note  that  any  eigen- 
value of  B which  has  a non-zero  eigenvalue  is  expressible  as 

e zie1. 


« 


A A 


f 


l 


I 


Let  <l>2,  ....  be  the  vectors  l z^  where  (zj, 

T ' 

zin)  is  in  the  null  space  of  AA  and  <}>  ^ * • ••»  ❖ n m 

be  the  vectors  orthogonal  to  the  in  direction  vectors  £ ^ 

...»  f,m-  Also  let  pj,  p2>  ...»  pn_m,  be  the  eigenvectors  of 
B (C)  with  an  eigenvalue  of  zero  (g2).  Then  it  follows 
U’j.  'i,2*  ....  <j>m_m.»  4'1 » ...»  4>n  m)  forms  a basis  for 

(P|,  p2»  ...»  P n _ in  * ) • It  should  be  noted  that  i|i^  is  orthogonal 
to  ij>.,  i = 1,  2,  ....  in-111'  and  j = 1,  2,  n-m'. 

J 

Now  we  consider  the  situation  where  two  sample  covariance 

matrices  C and  C are  generated  from  n and  n + 1 antenna  ele- 

ments  and  let  the  gramians  of  the  sets  of  direction  vectors 

be  A and  A respectively.  It  is  clear  that  A = A + }:*}: 
n n+ 1 n n + 1 n 

where  £ =(o^»  o!,'.  ...»  o(")  and  hence  that  the  eigenvectors  of 

C and  C are  in  general  different.  Because  of  the  nature 
of  the  direction  vectors  it  is  also  clear  that  the  two  vectors 
v.  and  v..  formed  from  the  first  and  last  n components  of  an 

i c.  2 i 

eigenvector  of  C|i  + j with  an  eigenvalue  greater  than  g belong 
to  the  linear  space  spanned  by  the  m direction  vectors  (with 
n components).  This  suggests  our  algorithm.  Generate  eigen- 
vectors for  C and  C and  select  those  eigenvectors  whose 
n n + 1 2 

eigenvalues  are  greater  than  T (>g  From  the  thus  selected 
set  of  eigenvectors  of  C[1  + 1»  generate  two  sets  of  vectors 
V i»  V2  by  choosing  the  first  and  last  n components  from  each 
eigenvector.  Augment  the  selected  eigenvectors  of  Cn  with 
vectors  from  and  V ^ using  Gram-Schmidt  orthogonal i zation 
procedure.  Once  this  is  done  we  apply  the  techniques  de- 

■ 

veloped  in  2.1  to  obtain  estimates  for  angular  locations. 
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2.3  Effect  of  Errors  on  Computational  Accuracy 


There  are  two  primary  sources  of  error  which  affect 
accuracy  of  estimates  of  angular  locations.  The  first  source 
is  the  A/D  converter  which  introduces  errors  due  to  finite 
levels  of  quantization.  The  second  source  is  the  antenna 
elements  whose  internal  noise  waveforms  may  not  be  uncor- 
related white  Gaussian  noises  with  identical  variances  after 
a finite  number  of  samples.  The  effect  of  this  error  is 
that  the  antenna  internal  noise  covariance  matrix  will  not 
be  a diagonal  matrix  with  identical  elements  but  a matrix 
with  unequal  diagonal  elements  and  non-zero  off-diagonal 
el ements . 

A simple  first  order  perturbation  analysis  (WILK  65) 
yields  useful  insight  into  the  effect  of  errors. 

Initially,  we  will  only  consider  the  effects  of  errors 
on  eigenvectors  and  show  how  to  minimize  them.  Let  A be 
a Hermitian  matrix  with  eigenvectors  x^,  x^»  . ...  xn  whose 
corresponding  eigenvalues  are  Xp  X2>  ....  Xn.  For  the 
perturbed  matrix  (A  + t:B)  let  x^e)  and  X^(e),  i = 1,  2, 

....  n be  the  eigenvectors  and  eigenvalues.  For  a first 
order  perturbation  analysis  we  assume  X^(e)  = X^  + k^e  and 
x.(e)  = x . + e { l sijxj)*  Then  it  follows: 
ij'j 

(A  + eB)xi (e)  = Xj (e)xi (e) 
collecting  first  order  terms  in  e we  have: 
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A{£  si jxj> + ^ - x.{  i sijXj}  + kiXi 
j i 


£ sij(xj  - Xt>xj  + Bx>  ■ kixi 
i T4  j 


Since  for  a Hermitian  matrix  we  can  choose  x^,  x^,  ....  xn 
to  be  orthonormal  we  have 

su  = xj8xi/(xi  - xj) 

i.e.,  xf(e)  = x(  * e I sjixj/(xi  - xj> 

j^i 

where  0..  = x?Bx..  Thus  the  eigenvectors  are  sensitive  to 
■ J £ J 

the  separation  of  eigenvalues. 

In  our  case  it  may  be  noted  that  we  are  interested  in 
eigenvectors  that  are  linear  combinations  of  direction  vectors 
and  concerned  only  when  the  vectors  orthogonal  to  the  direc- 
tion vectors  can  contaminate  the  desired  eigenvectors.  This 
immediately  suggests  that  in  general  the  larger  the  eigen- 
value, the  more  accurate  its  corresponding  eigenvector  will 
be.  Thus  in  the  algorithm  described  in  2.2  where  multiple 
eigenvector  sets  are  computed'  and  a threshold  T is  set  to 
select  eigenvectors,  T reflects  the  accuracy  of  eigenvectors 
containing  the  desired  information. 
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v wyi 


Now  we  consider  the  effects  of  errors  on  the  roots  of 
the  principal  polynomial  G(x)  using  first  order  perturbation 
analysis.  We  assume  that  the  obtained  eigenvectors  of  the 
combined  signal  covariance  matrix  C is  orthonormal  and  Coral- 
lary  2 is  applied  to  find  G(x).  We  indicate  the  first  order 
errors  in  the  eigenvectors  as 


and  in  X as  X + eXe  and  obtain  the  following  perturbed  princi- 
pal polynomial  Gp: 

Gp ( x ) = G(x)  + eGe(x) 

where  Gg(x)  = (w[eWf  + w£w£)  Vk_j(x)  + Xexk_1.  Note  the  co- 
efficients of  Gg(x)  cannot  exceed  in  magnitude  21.  Now  we 
consider  the  roots  of  Gp  ( x ) and  their  relation  to  those  of 
G ( x ) . We  have  by  Taylor's  theorem 

Gp ( x+A ) = G(x+A)  + eGe(x+A) 

= G(x)  + AG'(x)+  eGe(x) 

+ eAG^(x) 

If  xR  and  xR+A  are  roots  of  G(x)  and  Gp(x)  respectively,  we 
have 


W + eW 


"k  + ewke 
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2.4  Computational  Analysis  for  Multiple  Source  Location 


In  this  section  we  perform  a computational  analysis  of 
the  algorithms  presented  in  preceding  sections  in  order  to 
assess  the  feasibility  of  implementing  digital  multiple 
source  location  in  real  time.  These  algorithms  have  the 
following  three  phases  of  computation: 

1.  E i gen va 1 ue/ E i gen  vector  Computation 

2.  Polynomial  Formation 

3.  Polynomial  Solution  on  the  Unit  Circle 

For  computational  analysis  we  select  two  specific  algorithms 
which  are  described  in  the  following.  In  the  first  algorithm 
only  one  e i gen va 1 ue/e i gen vec tor  computation  is  performed  (for 
an  nxn  Hermitian  matrix)  and  the  principal  polynomial  G(x)  is 
computed  and  solved  on  the  unit  circle.  To  find  the  zeros  of 
the  principal  polynomial  on  the  unit  circle,  the  quantities 
| G [ exp  ( j k 0 ) ] | 2 , G = 2ir/KD,  k = 0,  1,  ....  k^-l,  are  computed 
to  locate  the  approximate  ranges  of  the  minima  of  | G ( x ) | 2 . 
Then  a Fibonaccian  scheme  with  KF  searches  (ORTE  70)  is  con- 
ducted  to  locate  minima  of  |G(x)|  and  hence  zeros  of  G(x). 

In  the  second  algorithm  which  is  more  general  than  the 

first,  eigenvalue/eigenvector  computations  are  performed  for 

nxn  and  (n  + 1)  x (n  + 1)  Hermitian  matrices.  Eigenvectors 

2 

whose  eigenvalue'  exceed  a preset  threshold  T(>g  ) are 
selected  and  Gram-Schmidt  orthogonal i zati on  procedure  Is 
performed.  Then  principal  and  auxiliary  polynomials  are 
computed  and  least  squared  error  solutions  are  determined 
( Section  2.1). 
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lor  o i gen va  1 ue/e  1 genvec tor  computations,  one  may  use 
Householder's  transformations  to  trl -dl agonal  1 ze  the 
Hermitian  matrix  and  apply  QR  transformations  to  obtain 
eigenvalues/eigenvectors  (W1LK  6b,  DAHL  / 4 ) . Assuming 
two  QR-i terations  per  eigenvalue  (DAHL  74)  an  eigenvalue/ 
eigenvector  computation  for  an  nxn  Hermitian  matrix  re- 
quires 20n^/3  + O(n^)  complex  multiplications  or  approxi- 

3 

mately  28n  real  multiplications.  For  polynomial  solution 
in  the  first  algorithm,  we  require  approximately  (KQ  + Kp)n 
complex  or  4(K^  + K(  )n  real  multiplications.  Assuming  that 
the  algorithm  has  to  be  performed  R times  every  second,  we 
require  a processing  rate  of  (28n  + 4(KQ  + Kp)n]R  multi- 

plications/second for  real  time  implementation.  (We  re- 
tain the  linear  term  in  n since  is  usually  large,  say 
around  1000).  Figure  2.4-1  shows  this  processing  rate  in 
MMPS  (million  multiplications/second)  against  the  number 
of  elements  for  • 1000,  K(r  - 60  and  R = 1,  10,  100f  1000. 

for  the  second  algorithm  a similar  computational  analysis 
reveals  that  we  require  approximately  [68n  + (K^  + Kj). 

(2  + kpriy)  n^JR  real  multiplications/second  if  the  algorithm 
has  to  be  performed  R times  every  second.  ^trig  rat^° 

of  computing  times  to  perform  the  sine  function  and  multiplica- 
tion. Figure  2.4-2  shows  the  required  processing  rate  against 
the  number  of  elements  for  + Kp  = 1060,  3 8 and  R - 

1,  10,  100,  1000. 

From  our  computational  analysis  and  the  present  and 
future  trends  of  computing  technology  (TURN  74),  It  is 
clear  that  real  time  implementation  of  digital  multiple 
source  location  procedures  is  feasible  within  this  decade. 

(The  only  technological  obstacle  may  lie  in  economic  con- 
struction of  A/D  converters  with  required  precision).  Also 
the  algorithms  underlying  the  procedures  are  amenable  to 


parallel  processing  and  Appendices  2 and  3 show  how  to  ex- 
ploit parallelism  inherent  in  matrix  operations  on  a par- 
ticular parallel  processing  computer  described  in  Appendix 
1 to  speed  up  computation.  A simple  computational  analysis 
shows  that  for  the  algorithms  proposed  for  multiple  source 
location,  a speed  up  proportional  to  p may  be  expected  using 
parallel  processing  where  p is  the  number  of  processors  in 
the  computing  system  and  is  in  the  range  20  - 100. 


3.  COMPUTER  IMPLEMENTATION  AND  NUMERICAL  EXPERIENCE 


A FORTRAN  program  was  written  to  implement  the  multiple 
source  problem  on  a PDP-11/35.  Only  the  first  algorithm  de- 
scribed in  2.4  has  been  completely  implemented  and  in  this 
section  we  report  on  the  numerical  experience  gained  by  using 
the  algorithm.  The  signal  sources/ jammers  and  the  internal 
antenna  noises  are  simulated  by  generating  independent  normal- 
ly distributed  random  numbers  (using  Box-Muller's  transforma- 
tion on  uniformly  distributed  numbers  in  [0,  1]  (DAHL  74)) 
with  zero  mean  and  prespecified  variances.  The  program  was 
written  in  double  precision  arithmetic  (56  bits  of  mantissa, 

8 bits  of  exponent)  and  can  handle  up  to  16  antenna  elements 
and  10  sources.  A listing  of  the  program  is  attached  as 
Appendix  4 of  this  report. 


Various  numerical  experiments  were  conducted  to  observe 
the  behavior  of  the  algorithm  in  terms  of  its  accuracy  in  de- 
termining signal  source  locations.  In  the  first  experiment, 
a linear  array  of  4 antenna  elements  with  interelement  spac- 
ing of  half  wavelengths  is  assumed  and  two  sources  are  simulated 
at  angular  locations  .25it  and  -.3n  radians.  The  signal  source 
variances  are  assumed  to  be  10  power  units  and  the  internal 
antenna  noise  power  in  each  element  is  taken  to  be  1 unit. 

Table  3-1  shows  the  accuracies  obtained  as  the  number  of 
samples  is  varied  from  10  to  100.  The  inaccuracies  in  esti- 
mates of  angular  locations  are  due  to  the  fact  that  the  in- 
ternal antenna  noise  covariance  matrix  is  not  an  identity 
matrix  as  assumed.  (In  the  simulation  if  the  internal  noise 
covariance  matrix  is  assumed  to  be  an  identity  matrix,  exact 
angular  locations  were  obtained.)  Table  3-2  shows  the  internal 
antenna  noise  covariance  matrices  after  10  and  50  samples  and 
the  eigenvalues  of  these  matrices. 
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Table  3-1.  Experimental  results  illustrating  the  effects  of  number  of 
samples  on  angular  accuracy. 


Table  3-2.  Antenna  internal  noise  covariance  matrices  and  their  eigenvalues 
after  10  and  30  samples. 


A * 


In  the  next  experiment,  the  effect  of  the  number  of 
antenna  elements  on  the  accuracy  of  estimates  is  investi- 
gated. The  number  of  samples  is  fixed  at  30  and  three 
signal  sources  are  simulated  at  . lir,  .2 u and  . 3 ir  radians 
with  mean  power  levels  at  10  units.  The  internal  antenna 
noise  power  is  again  fixed  at  1 unit.  Table  3-3  shows 
the  results  obtained  as  the  number  of  elements  is  varied 
from  4 to  8.  One  can  immediately  see  that  increasing  the 
number  of  antenna  elements  results  in  improved  estimates 
of  angular  locations. 

In  the  final  experiment  the  angular  resolution  capa- 
bilities of  the  algorithm  are  observed  as  the  source  sig- 
nal powers  and  the  number  of  antenna  elements  are  changed. 

Two  signal  sources  are  assumed  and  the  number  of  samples 
is  fixed  at  30.  Initially  the  signal  powers  are  fixed 
at  10  units  and  the  angular  locations  at  .In  and  ,2tt 
radians.  Table  3-4  shows  the  results  of  this  experiment. 

As  the  second  source  is  moved  closer  to  the  first,  the 
algorithm  coalesces  both  sources  into  one.  (The  extrane- 

p 

ous  minima  of  |G(x)|  on  the  unit  circle  are  marked  by 
asterisks).  Thus  one  can  see  from  the  table  that  when 
the  sources  are  located  at  .lit  and  .I3tt  or  closer,  the 
algorithm  cannot  resolve  the  sources  satisfactorily. 

However,  increasing  the  power  levels  of  the  sources  or 
the  number  of  elements  results  in  the  sources  being  re- 
solved. 

In  Table  3-5  we  show  an  example  of  coalescing  two 
sources  into  one  and  the  results  of  choosing  various 
thresholds.  Two  sources  are  simulated  at  .1  and  . 1 2 it 
with  power  levels  fixed  at  10  units.  As  before  the  in- 
ternal antenna  noise  power  is  assumed  to  be  1 unit  and  a 
four  element  linear  array  is  assumed.  By  setting  threshold 

<:  i 
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antenna  elements  on  angular  accuracy. 


Antenna  Array  Interelement  Spacing:  Half  Wavelength 
Number  of  Signal  Sources:  2 
Internal  Antenna  Noise  Variance:  1 
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at  various  levels,  the  number  of  eigenvectors  selected  is 
varied  from  one  to  three.  The  extraneous  minimas  are  marked 
with  an  asterisk.  The  same  experiment  is  repeated  with 
sources  located  at  . 2 5tt  and  -.3tt  and  as  can  be  noted,  no 
extraneous  minima  appeared  when  three  eigenvectors  are 
selected.  Though  extraneous  minima  can  be  eliminated 
theoretically  by  considering  auxiliary  polynomials,  it  is 
not  clear  whether  such  is  the  case  when  the  data  is  con- 
taminated by  error  noise. 

In  order  to  test  the  effectiveness  of  the  algorithm 
when  data  was  obtained  from  an  actual  antenna  array,  mag- 
netic tapes  were  obtained  from  NRL  through  the  courtesy  of 
Mr.  Fred  Staudaher  and  processed.  These  tapes  contained 
covariance  matrices  obtained  from  a 4-element  linear  array. 
The  number  of  sources  varied  from  0 to  10  and  an  analog-to- 
digital  conversion  scheme  with  512  samples  was  used.  To 
compute  the  covariance  matrices,  1024  samples  were  used. 

The  covariance  matrices  for  the  cases  where  the  number 
of  sources  is  from  1 to  3 were  processed  and  the  results 
are  shown  in  Tables  3-6  to  3-9.  As  can  be  seen,  the  re- 
sults are  encouraging  for  one  source  where  the  angular 
errors  are  below  4 degrees  (Table  3-6).  For  two  sources 
except  for  two  cases  (File  Numbers  17  and  18  in  Table  3-7), 
the  angular  errors  are  below  5 degrees.  For  File  Numbers 
17  and  18,  it  appears  that  the  algorithm  coalesces  two 
sources  into  one.  These  cases  and  File  Numbers  14  and  15 
were  run  again  selecting  different  numbers  of  eigenvectors 
for  polynomial  formation  (Table  3-8).  (The  extraneous 
minima  as  before  are  marked  with  an  asterisk).  For  the 
three  source  cases  (Table  3-9),  the  algorithm  appears  to 
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Number  of  Antenna  Elements: 


Table  3-6.  Results  of  processing  NRL  magnetic  tapes  wi 


.-.tier  of  Ar. tenna  Elements: 


able  3-7.  Results  of  processing  NRL  magnetic  tapes  with  two  sources 


Nuaoer  of  Antenna  Elements: 
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Table  3-9.  Results  of  processing  NRL  magnetic  tapes  with  three  sources. 


perform  less  satisfactorily.  These  cases  were  run  under 
the  assumption  of  one,  two  and  three  sources  and  as  Table 
3-9  shows,  in  only  one  case  (File  Number  22)  were  the 
three  sources  resolved  satisfactorily. 

The  main  source  of  computational  errors  is  in  the  A/D 
converter  which  uses  only  512  levels.  Hence  the  data  is 
accurate  only  up  to  the  second  decimal  digit  and  this  seems 
to  introduce  significant  errors  in  estimates.  Also,  the 
number  of  antenna  elements  (4)  is  quite  small  to  resolve 
sources  close  to  each  other  satisfactorily.  It  is  hoped 
that  we  may  obtain  in  the  near  future  data  from  an  antenna 
receiver  system  with  more  elements  and  precision. 
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4.  CONCLUSIONS  AND  SUGGESTIONS  FOR  FUTURE  RESEARCH 

This  report  presented  algorithms  for  locating  simultaneous- 
ly multiple  signal  sources  by  means  of  linear  arrays  and  digital 
signal  processing.  Computational  complexity  analysis  of  these 
algorithms  shows  that  it  is  possible  to  implement  them  in  real 
time  on  high  speed  parallel  processing  computer  systems  with- 
in the  next  decade.  Simulations  of  the  algorithms  yield  satis- 
factory results  after  a small  number  of  input  samples  (30-50) 
and  preliminary  results  indicate  that  one  can  be  optimistic  in 
expecting  high  resolution  and  good  accuracy  when  a relatively 
small  number  (16-32)  of  antenna  elements  is  used. 

At  this  time  it  is  appropriate  to  compare  the  proposed 
algorithms  with  that  of  Berni  (BERN  75).  For  Berni's  algorithm, 
one  is  forced  to  restrict  the  number  of  elements  to  m+1  where 
m is  the  number  of  sources  to  be  located.  As  indicated  in 
Section  3,  better  accuracies  and  resolution  can  be  expected 
by  having  a large  number  of  antenna  elements  (refer  to  Tables 
3-3  and  3-4)  and  thus  Berni's  algorithm  may  not  be  preferable 
in  this  sense.  Also  it  may  be  clear  from  the  description  of 
his  algorithm  that  it  is  not  applicable  when  the  sources  are 
correlated  and  the  signal  source  covariance  matrix  is  singular. 

For  future  research,  it  is  suggested  that  a thorough 
experimental  investigation  be  conducted  to  study  the  behavior 
of  the  proposed  algorithms  in  terms  of  accuracy,  resolution 
and  convergence  at  various  signal  source  power  levels  and 
angular  locations  and  at  various  numbers  of  signal  sources. 

The  future  research  should  concentrate  on  evaluating  the 
effects  of  errors  (due  to  A/D  conversion  as  well  as  the  non- 
whiteness of  antenna  internal  noise)  on  computational  accuracy 
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and  methods  of  minimizing  these  effects.  The  intriguing 
problem  of  extraneous  minima  (mentioned  in  Section  3)  should 
also  be  examined  to  see  whether  they  can  be  completely  elimin 
ated.  The  effectiveness  of  the  proposed  algorithms  in  locat- 
ing signal  sources  that  have  a singular  covariance  matrix  is 
an  area  that  is  yet  to  be  explored  and  the  future  research 
may  focus  on  the  methods  of  using  the  proposed  algorithms  to 
combat  the  tracking  errors  caused  by  multipath  propagation. 

It  would  also  be  of  interest  to  investigate  whether  circular 
and  conformal  arrays  may  be  used  instead  of  linear  arrays  for 
multiple  signal  location.  Finally  the  problem  of  recovering 
multiple  signals  with  different  angles  of  arrival  should  be 
of  major  interest  in  the  fields  of  communication  and  defense 
and  hence  be  pursued  in  the  future. 


t 


■ 


39 


REFERENCES 


BERN  75 

Berni  , A.J.,  "Target  Identification  by  Natural  Resonance 
Estimation",  IEEE  Trans.  Aerospace  and  Electronic  Systems, 
Vol.  AES-11,  No.  2,  March  1975,  p.  147-154. 

DAHL  74 

Dahlquist,  6.,  et  al . , Numerical  Methods,  Prentice-Hall, 

NY,  1974. 

GANT  60 

Gantmacher,  F.R.,  The  Theory  of  Matrices,  Vol.  1,  Chelsea, 
NY,  1960. 

NERI  67 

Nering,  E.D.,  Linear  Algebra  and  Matrix  Theory,  Wiley,  NY, 
1967. 

WILK  65 

Wilkinson,  J.H.,  The  Algebraic  Eigenvalue  Problem,  Clarendon 
Press,  Oxford,  1965. 


' " T""  1 


Glossary  of  Symbols 

No.  of  Antennas:  n 

No.  of  Signal  Sources:  m 

Signal  Sources/ Jammers : J,_  J_,  ...,  J 

1 d m 

Angular  Locations  of  Sources:  Qj.  02»  •••*  9m 

Wavelength  of  the  Carrier  Signal:  v 

Direction  Vectors:  E,  , £„ C 

I d m 


1 

exp  { j ( 2-rrd/ v ) sin  0^} 
exp  { j (4-rrd/v)  sin  0^ } 

I 

I 


exp  { j(n-l)  (2ird/v ) sin  0.} 


“ ~ ^ 1 1 ^2  I * * * 


Complex  Envelope  Amplitudes  of  Signal  Sources: 
Sl*  S2 * sm 

Source  Signal  Covariance  Matrix: 

A - I I ai j | | - | ! E(s*Sj ) 1 1 
Received  Signal  Covariance  Matrix: 

8 ' 1 1 bi  j ! I ' I Wl 

i »j 

Antenna  Internal  Noise  Covariance  Matrix: 
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Combined  Signal  Covariance  Matrix: 
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= g i 

+ B 

Eigenvectors  of  C 

: 

» C 2 » 

....  ci 

P1 

» p£  » 

• • • » P j 

Eigenvalues  of  B: 

*1- 

3 2 * 

• • » 3 ^ » 

Eigenvalues  of  C: 

V 

y2»  ■ 

...  yv 
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KEY  FEATURES  OF  THE  G-471,  A NEW  PARALLEL/ASSOCIATIVE  COMPUTER 


1. 


INTRODUCTION 


The  G-471  is  a 4th  generation  parallel/associative 
computer  whose  architecture  and  technology  produce  signifi- 
cant improvements  in  performance  and  maintainability  over 
earlier  designs,  at  considerably  lower  cost.  It  is  des- 
cribed in  more  detail  in  Chapter  2. 

Its  most  important  advantages  are: 

Extremely  high  processing  rates  (64  - 4096  MIPS) 

Extremely  high  memory  bandwidth  (512  Mbytes/sec  for 
every  32  Processing  Elements,  16,384  Mbytes/stc  for 
1024  PEs) 

Main-stream  technology 

Much  lower  cost  than  conventional  scientific  computers 
Minimal  learning  time  for  programming 
Microprogrammabi 1 i ty 

Higher-order  languages  (FORTRAN,  BASIC,  COBOL) 
Multiprocessing  for  higher  speed  and  resource  utilization 
Parallel  pipeline  capability 

Operation  as  a stand-alone  unit  or  as  a peripheral  to 
other  computers 

Operation  with  off-the-shelf  disks  and  tapes  as  mass 
memory 

High  reliability  through  fault  tolerance  and  easy 
maintainability. 
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2.  DETAILED  DESCRIPTION 

The  figure  on  the  following  page  shows  the  block 
diagram  of  the  G-471.  The  Sequential  and  Control  Computer 
(SCC)  is  a standard  or  militarized  minicomputer  or  a larger 
machine  if  desired.  It  controls  the  operations  of  the  Proc- 
essing-Element and  Data-Routing  Element  Arrays  (PE&DREA),  and 
it  also  performs  sequential  computations.  The  optimum  size 
of  this  Sequential  and  Control  Computer  depends  on  the  size 
of  the  array  to  be  controlled  (16  to  1024  Processing  Elements); 
on  the  number  of  control  actions  required  (if  the  processing 
elements  perform  mostly  routine  operations  which  are  programmed 
in  the  PE  local  memories  the  demand  on  the  control  computer 
is  small;  if  the  PE  programs  require  frequent  changes,  the 
demand  on  the  control  computer  is  higher);  and  on  how  many 
sequential  operations  must  be  performed  in  addition  to  the 
array-control  functions.  The  Sequential  and  Control  Computer 
can  be  supplied  with  the  overall  system,  or  it  can  be  a cus- 
tomer-owned computer  for  which  only  the  necessary  interfaces 
to  the  PE  and  Data-Routing  Element  Arrays  and  the  mass-memory 
controllers  are  supplied.  The  Sequential  and  Control  Computer 
has  standard  peripherals  such  as  a CRT  terminal  and  high-speed 
printer,  and  it  normally  runs  under  a real-time  disk-operating 
system . 

The  Processing  Elements  (PEs)  in  the  PE  Array  contain 
standard  minicomputers,  each  of  which  processes  16  bits  in 
parallel.  The  local  memory  associated  with  each  PE  is  expand- 
able to  at  least  56  kbytes.  This  storage  area  can  be  assigned 
to  data  or  program,  RAM  or  ROM,  in  any  mix. 

The  PE  designs  have  been  enhanced  significantly  as  com- 
pared to  standard  minicomputers  to  make  the  most  frequent 
operations  in  the  PE  Array  particularly  fast.  For  example, 

DMA  (direct-memory  access)  to  portions  of  the  PE  memory  is  pos- 
s i 6 1 e while  other  portions  are  simul taneously  accessed  by  the  CPU. 
Each  microprogrammabl e high-speed  arithmetic  enhancement  pro- 
cessor carries  out  32-bit  floating-point  multiplications  at  a 
rate  of  4 MIPS  (128  MIPS  for  32  PEs). 

PEs  can  be  readily  paralleled,  with  arrays  ranging  from 
16  to  1024  PEs  (256  to  16,384  bits  in  parallel). 

Each  PE  can  directly  address  up  to  16  Mbytes  of  semiconductor 
RAM  Central  Working  Storage  (CWS).  The  CWS  is  partitioned  into 
at  least  as  many  memory  banks  as  there  are  PEs  to  permit  parallel 
access.  Each  PE  and  each  Memory  Bank  is  provided  with  a 64-bit 
wide  memory  transfer  channel  allowing  a transfer  rate  per  channel 
of  16  Mbytes/sec.  Each  group  of  32  PEs  therefore  has  a memory 
bandwidth  of  512  Mbytes/sec. 
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Each  PE  also  has  access  to  the  mass  memory,  typically  a 
bank  of  disks,  using  controllers  which  are  operated  by  the 
Sequential  and  Control  Computer. 

In  some  applications  the  data  to  be  processed  may  not  be 
supplied  from  mass  memory  but  from  real-time  sensors  with  a 
high  data  rate,  such  as  video  sensors,  radar,  sonar,  large 
distributed  sensor  arrays  etc.  For  this  purpose  Real-Time  I/O 
Channels,  each  with  a transfer  rate  of  16  Mbytes/sec,  are  pro- 
vided . 

The  Data-Routi ng-El ement  Array  performs  the  communication 
function  among  the  PEs,  the  Real-Time  I/O  Channels,  the  CWS 
Memory  Banks  and  the  Mass-Memory  Modules.  It  is  basically  a 
programmable  cross-point  switch  whose  switch  settings  are  deter- 
mined dynamically  by  the  memory  addresses  from  the  PEs,  the  Real- 
Time  I/O  Channels,  the  Mass-Memory  Controllers,  and  the  Control 
Computer.  The  DRE  Array  operation  is  transparent  to  the  programs. 

The  overall  operation  of  the  PE  and  Data-Routi ng-El ement 
Arrays  is  controlled  by  the  Control  Computer.  The  control 
functions  are  as  follows: 

a.  Activate  any  single  PE  or  any  group  of  PEs. 

b.  Transmit  instructions  to  any  PE,  singly  or  in  a 
broadcast  mode  (instructions  can  also  be  loaded 
from  CWS  and  Mass  Memory). 

c.  Transmit  data  to  the  PEs  (in  the  majority  of  appli- 
cations most  of  the  data  would  come  from  the  Mass 
Memory  and  from  the  Real-Time  I/O  Channels). 

d.  Receive  data  and  "response  signals"  from  the  PEs. 

e.  Control  the  system  reconfiguration  to  achieve  fault 
tolerance . 

f.  Allocate  tasks  and  subtasks  to  available  processing 
el ements . 

g.  Set  up  protection  registers  to  limit  the  PEs  to 

permissible  portions  of  the  Central  Working  Storage  (CWS). 

A striking  feature  of  the  design  is  the  large  amount  of 
off-the-shelf  hardware  used,  specifically  the  control  computer, 
the  disk  controllers,  the  mass  memory  (disk,  drum,  tape),  the 
control  console  and  printer,  and  the  processing  elements  them- 
selves (except  for  the  high-speed  hardware  enhancement  options). 

A major  side  benefit  of  this  architecture  is  that  it  runs 
under  standard  real-time  disk  operating  systems  and  has  all 
higher  language  capabilities,  as  well  as  console  and  printer 
software.  Remote  access  via  modem  could  also  be  standard. 
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3.  ADVANTAGES  OVER  EARLIER  DESIGNS 

Extremely  High  Processing  Rates  (64  - 4096  MIPS) 

The  minimum  configuration  of  16  PEs  will  process  up  to 
256  Mbytes/sec.  The  corresponding  number  for  the  "jumbo" 
1024-PE  array  is  16384  Mbytes/sec. 

Extremely  High  Memory  Bandwidth  (512  Mbytes/sec 

Per  32-PE  Array*] 

The  G-471  design  recognizes  that  the  high  processing  power 
of  the  PEs  must  be  supported  by  the  ability  to  move  data  very 
rapidly  in  and  out  of  a large  Central  Working  Storage  (CWS). 
Thus  it  is  possible  for  each  PE  to  directly  address  16  Mbytes 
of  CWS  and  to  transfer  data  at  a rate  of  16  Mbytes/sec.  Since 
PEs  can  operate  in  parallel,  the  effective  memory  bandwidth 
for  a 32-PE  Array  is  512  Mbytes/sec. 

Main  Stream  Technology 

The  G-471  is  the  first  parallel/associative  processor  in 
main-stream  technology  using  off-the-shelf  microprocessor 
assemblies,  minicomputers,  controllers  and  disks.  The  amount 
of  custom  circuitry  is  minimized.  The  benefits  of  millions  of 
dollars  of  minicomputer  R&D  in  hardware,  software,  diagnostics 
and  education  are  incorporated  into  the  G-471. 

Low  Cost 


Because  it  utilizes  many  fully  designed  and  debugged 
system  modules  which  are  in  mass  production,  its  hardware  and 
software  development  costs  are  drastically  lower  than  a machine 
designed  from  the  gate  level  up.  With  the  packaging  schemes 
for  most  of  the  components  worked  out,  and  with  spares  readily 
available,  the  manufacturing  and  maintenance  costs  are  also 
greatly  reduced. 

The  cost  reduction  over  other  machines  with  similar 
capability  amounts  to  a factor  between  6 and  10. 

Minimal  Learning  Time  For  Programming 

The  G-471  is  an  associative/parallel  processor  computer 
which  "you  can  program  on  Day  One".  Rather  than  requiring  its 
own  fundamental  machine  language,  the  instruction  repertoire 
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of  the  6-471  is  derived  from  the  instruction  repertoires  of 
the  minicomputers  used.  Thus,  a programmer  familiar  with 
minicomputer  software  can  use  the  G-471  with  a minimum  of 
learning  time.  The  execution  times  of  all  instructions  are 
known . 

Microprogrammabil i ty 

The  processing  elements  in  the  array  are  individually 
microprogrammabl e , thus  allowing  the  ready  tailoring  of  the 
G-471  to  a given  application. 

Higher-Order  Languages 

Programs  written  in  FORTRAN,  BASIC,  and  COBOL  can  be 
run  on  the  G-471. 

Multiprocessing  for  Higher  Speed  and 

Resource  Utilization 


Different  processors  in  the  G-471  can  carry  out  different 
programs  simultaneously,  a major  difference  from  conventional 
parallel  processors  which  usually  perform  only  one  instruction 
at  a time  on  all  data  in  parallel. 

Fault  Tolerance  and  Easy  Maintainability 

The  special  architecture  of  the  G-471  with  its  high  degree 
of  modularity,  and  the  large  amount  of  diagnostics  available 
from  the  minicomputer  systems  components,  makes  the  G-471 
highly  fault  tolerant  and  maintainable,  producing  a much  great- 
er availability  than  is  customary  in  machines  of  this  complex- 
ity. 


Easy  Expandability 

The  G-471  is  expandable  in  four  directions: 

a.  Expansion  in  the  number  of  parallel  processing 
elements  (up  to  1024  PEs,  processing  1 6 K bits 
in  parallel). 

b.  Expansion  of  the  fast  local  semiconductor  memory 
associated  with  each  parallel  processing  element 
up  to  56  kbytes  (58  Mbytes  in  1024  PE  machine). 

c.  Expansion  of  Central  Working  Storage  to  16  Mbytes 
per  32  PEs. 
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d.  Expansion  of  Mass  Memory  by  adding  parallel  disks 
and  controllers. 

Operation  As  a S ta  n_d -A!  o n e Unit  Or  Asa  P e r i p her  a 1 

to  Other  Computers 

The  G-471  is  configured  to  operate  as  a stand-alone 
unit,  or  the  Process i ng- El emen t and  Data-Routing-Element 
Array,  the  Central  Working  Storage  and  the  Mass  Memory  can  be 
interfaced  with  another  computer  to  operate  as  a peripheral 
which  carries  out  specific  parallel  or  associative  functions. 

The  software  is  separable  in  the  same  way  as  the  hardware.  This 
is  of  advantage  when  the  associative  capability  is  to  be  added 
to  an  operational  computer  system  with  a large  amount  of  other 
fully  debugged  software  which  must  not  be  disturbed. 

Operation  With  Off-The-Shelf  Disks  As  Mass  Memory 

The  G-471  has  been  configured  to  operate  with  standard 
off-the-shelf  disk  units  as  mass  memory.  This  has  the  advan- 
tage that  standard  disk  controllers  can  be  used  with  their  de- 
bugged diagnostics  and  maintenance  procedures.  It  also  allows 
programming  of  the  G-471  on  other  machines,  and  the  transfer 
of  the  programs  via  removable  disk  cartridge. 
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4.  TYPICAL  APPLICATION  AREAS 

The  advantages  of  the  G-471  come  into  play  in  applications 
which  require  very  high  processing  rates  and  large  fast  memory, 
and  where  the  processing  algorithms  allow  a high  degree  of 
parallelism. 

A few  typical  examples  in  this  area  are: 

Image  (picture)  processing,  especially  restoration 
and  enhancement,  including  non-standard  algorithms. 

Radar  signal  processing. 

Sonar  signal  processing. 

Fast  Fourier  Transforms  (FFT)  and  other  array  processes. 
Real-time  beam  forming. 

Post  processing  of  image,  radar,  sonar  signal -processor 
output  data. 

Seismic  signal  processing. 

Correlation  of  signals  from  large  distributed  sensor 
arrays . 

Weather  prediction. 

EW/ECM  problems  which  involve  high  data  rates. 

Mapping  algorithms. 

Character  string  (associative)  searches  of  large  data 
bases . 

Multicriteria  searches  of  large  files. 

Mi ss i 1 e- tracki ng  and  air  traffic  control  problems. 
Monte-Carlo  analysis. 

High-speed  simulation. 

Partial  differential  equations. 

Matrix  manipulations. 

High-speed  simulation. 

The  G-471  is  not  only  a data  processor  but  can  also  be 
used  as  a large  array  controller  since  the  outputs  from  the 
individual  PEs  can  individually  control  a variety  of  external 
devices. 
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5.  RELIABILITY/MAINTAINABILITY 


The  G-471  has  been  specifically  designed  for  high  reliability, 
maintainability  and  availability,  using  a combination  of  per- 
formance monitoring,  fault  location,  self-repair,  redundancy 
and  on-line  maintenance.  The  design  is  taking  advantage  of  many 
years  of  pioneering  expertise  in  fault-tolerant  digital  system 
development  at  W.  W.  Gaertner  Research,  Inc. 

The  Control -Computer  Complex  can  be  made  dual  or  triple 
redundant,  with  voting  if  necessary.  The  mass  memory  (bank 
of  disks)  and  its  controllers  are  electrically  and  mechanically 
segmented  allowing  for  redundant  storage  of  data  and  programs 
where  desirable.  The  PEs  and  their  memories  as  well  as  the  portion 
of  the  Data-Routing  Network  which  services  the  particular  PE  are 
mounted  on  individual  boards  which  can  be  electrically  bypassed 
and  mechanically  removed.  The  Central  Working  Storage  is  par- 
titioned into  separate  memory  banks  and  it  can  be  electronically 
reconfigured  to  neutralize  hardware  faults. 

Performance  monitoring  and  fault  location  software  are 
provided,  and  on-line  maintenance  is  possible  (i.e.  replacement 
of  a faulty  module  can  take  place  with  the  rest  of  the  system 
operating) . 


6.  SIZING  OF  A SPECIFIC  G-471  SYSTEM 

A G-471  system  is  sized  by  analyzing  the  algorithms 
which  are  used  in  the  particular  application  area  (e.g.  image 
processing)  and  by  determining  the  required  maximum  execution 
rate  and  therefore  the  number  of  processing  elements,  the 
amount  of  high-speed  local  storage,  semiconductor  Central 
Working  Storage  and  Mass  Memory  (disk,  tape),  real-time  I/O 
facilities,  as  well  as  the  desired  peripherals  for  the  control 
computer. 

W.  W.  Gaertner  Research,  Inc.  is  available  on  a con- 
sulting basis  to  analyze  customer  applications  and  to  rec- 
ommend an  optimum  configuration  for  a G-471. 
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1.  MATRIX  MULTIPLICATIONS 

i 

Matrix  multiplications  possess  a large  degree  of  paral- 
lelism and  lend  themselves  to  parallel  processing.  Consider 
the  multiplication  of  two  square  matrices  A and  B of  order 
nxn.  Partition  the  matrices  into  k submatrices: 


" 
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A1B2 
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A2 

A2B  1 

A2B2 

-- 
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L 1 1 1 -1 
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• 0 
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O 

0 

0 

• 

• 

o 

0 

• 

< 

1 

AkBl 

AkB2 

-- 

AkBk 

2 

It  can  be  seen  that  k independent  matrix  multiplica- 
tions have  to  be  performed  and  each  matrix  multiplication 
3 2 

requires  n /k  multiplications  and  additions. 

Assuming  there  is  no  conflict  in  accessing  A.  and  B., 

* J 

the  time  required  to  perform  the  matrix  multiplication  on 
p processors  is 


(n3/k2)  |~ k 2 / p~ J = n3/p  mtu 

where  mtu  (multiplication  time  unit)  is  the  time  required  to 
perform  a single  real  multiplication.  (We  estimate  the  com- 
putational complexities  in  terms  of  mtu's). 
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-> 


In  matrix  multiplication  the  main  problem  lies  in 
storage  organization  of  matrix  elements  to  ensure  con- 
flict-free access  rather  than  the  exploitation  of  com- 
putational parallelism.  In  the  following  we  consider 
multiplication  of  1000x1000  matrices  and  10,000x10,000 
matrices  to  illustrate  the  organizational  problems  one 
encounters.  To  simplify  the  presentation  of  algorithms 
we  assume  the  following  computer  configuration: 


CWS:  10  banks  B K . 


BK, 


BK10  of  256 


1 » u 2 * • • ' 

Kbytes  each.  (2  Mbytes  of  Central 

Storage ) . 

PEs-  10  PEs  each  with  a local  storage  of  48 
Kbytes  in  three  memory  banks. 

Disk  Units:  1-2  units  with  a transfer  rate 
of  10  Mbytes/second. 

1000  x 1000  Matrix  Multiplication 


Partition  the  matrices  as  follows 


1000  100  100  100 


2 

'A1  " 

'A181 

Al82 

-- 

aiBio" 

2 

A2 

A2B1 

A2B2 

A 2 B 1 0 

[B1  |B2 

— |Bio]  1000  ’ 

• 

• 

• 

• 

• 

• • 

O 

• 

• 

o 

2 

A500 

A10B1 

A 1 0 B 2 

-- 

A 1 08  1 0 
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1.  Load  B^,  B^,  ...  Bp  into  the  ten  banks  of  CWS. 

I 

2 . Set  i = 0. 

3.  Load  A.  + p App  •••*  A.  + p into  the  local  storage 
(PELS)  of  the  ten  processors. 


4.  Perform  10  computational  steps  described  in  the  follow- 
ing. (The  procedure  is  described  for  i = 0). 


P1 

P2 

P3 

P10 

step  1 

A1B1 

A2B2 

A3B3 

— 

A 1 0B  10 

step  2 

A 1 B 10 

A2B1 

A3B2 
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A 1 0 B 9 

• 

o 

• 

0 
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O 

o 

• 
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• 

O 

• 

0 

9 

• j 

o 

o 

6 

O 

• 

step  10 

A1B2 

A2B3 

A3B4 

— 

A10B1 

In  step  1,  P.  performs  A . B . and  stores  the  result  in 


AiBi-l 


the  bank  containing  B.;  in  step  2,  P . performs  ^mQd  p^ 

and  stores  the  result  in  the  bank  containing  B._j  ^ilQd  pj 
and  so  on. 


5.  Transfer  the  partial  products  accumulated  in  CWS  to 
the  disk  unit(s). 

6.  If  i = 490,  stop.  Otherwise  i *-  i + 10  ano  go  to  step  2. 
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Estimation  of  Time  Requirements 
# 

t^:  Time  in  nanoseconds  to  perform  a single  real 

multiplication  250). 

tc:  Time  in  nanoseconds  to  store  (fetch)  in  ( from) 

CWS  without  conflicts  (/^500). 

rfc:  Disk  transfer  rate  in  Mbytes/second  (/^10). 

It  is  assumed  2t^tj.  and  that  each  matrix  element  is 
8 bytes  long. 


To  perform  A.B.,  we  fetch  (1000x100)  matrix  elements 
■ 

and  store  (2x100)  elements.  We  execute  (2x1000x100)  multi- 
plications. Since  2 1 t c , the  process  is  computational 
bound . 


Time  "** MATMUL  recluired  f°r  the  entire  multiplication 
can  be  given  as 

^MATMUL  = + (24/rt)  seconds. 

Hence  for  typical  tM(A-»250)  and  rt(/v^l0)  we  have 

TMATMUL  ^ 30  seconds- 


In  general  if  we  have  p processors  and  equal  number  of  banks 
in  CWS, 


TMATMUL  = (tM/p)  + ^24/rt)  seconds . 


■ ; 


i 
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I 

3.  Perform  steps  a)  and  b)  in  parallel. 


a)  Ass  i gn*  processor  PE.  to  BK^  and  compute  A^.B^j 
(of  dimensions  x X y)  using  the  outer  product 
expansion: 


1000 


XY  = 


X1Y1  + X2Y2  + 


+ X1000  Y 1000 


b)  Transfer  A. # i+1  (mod  1Q)  from  the  disk  to 
BKi  + l (mod  10)  • 

Figure  1-1  shows  the  contents  and  assignment  of  the 
banks . 

4.  Perform  steps  a)  and  b)  in  parallel. 

a)  Assign  P.  tQ  BKi  + 1 (mQd  1Qj  to  compute  A^  i+J  (mod  1Q). 
Bi+1  (mod  10),  1 usin9  outer  product  expansion. 

b)  Transfer  A(>  ,+2  (n|od  10)  from  the  disk  to  BK,t2  (mod  10). 

Figure  1-2  shows  the  contents  and  assignment  of  the 
banks . 
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After  the  processors  compute  the  x X y products 
(shown  hatched)  proceed  to  the  next  lOx  rows  of  A and  re- 
peat the  process. 

* 

y 


After  exhausting  the  rows,  repeat  the  process  by 
selecting  next  y columns  until  B is  completely  processed. 

Determination  of  x and  y and  Time  Estimates 

x and  y are  determined  as  follows: 


Number  of  multiplications 
performed  by  a processor 
in  each  stage 


1000  xy 
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: 1000  px 


Number  of  memory  accesses 

> : 1000  (x  + y) 

necessary  to  compute  the  product  I 

Number  of  fetches  from  the  ] 

\ : 1000  px 

disc  to  buffer  the  elements  of  A I 

For  complete  overlap  of  disc  transfers: 

(1000  xy)  t^lOOO  (x+y)  tc  + [(1000  px)  t Q/ D ] . 

Here  D is  the  number  of  disks  and  tQ  is  the  time  in  nanoseconds 
to  transfer  a byte.  To  accommodate  three  submatrices  in  each 
bank : 


8(2x+y)  _>  Mg 

where  Mg  is  the  capacity  of  each  bank  in  Kbytes  (assumed  to 
be  256). 


For  tM  = 250,  t^  = 500,  tQ  = 800  and  10  processors,  com 
plete  overlap  of  disk  transfers  can  be  achieved  for  D = 3. 
For  this  case  x and  y can  be  given  as: 


21, 

22 

19, 

20 

18 
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and  time  for  matrix  multiplication: 


. ( 10,000) 
MATMUL 

3 

10  t^/p  seconds 

~ 

7 h o u r s . 

For  p = 20  and 

D = 6, 

.(10,000) 

MATMUL 

3.5  hours. 

For  D/p  < . 3 

(p:  no. 

of  processors) 

.(10,000) 

MATMUL 

= MAX 

/l03tM/p,  103(  tc  + t()/D)/y' 

For  x = 6,  y = 

20  and 

D = 1: 

.(10,000) 

MATMUL 

= 16 

hours . 
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APPENDIX  3 


MATRIX  EIGENVALUE/EIGENVECTOR  COMPUTATION 
ON  A PARALLEL  PROCESSING  COMPUTER 


•I 


3.  MATRIX  EIGENVALUE/EIGENVECTOR  COMPUTATION 

We  consider  real  symmetric  matrices  and  analyze  the 
speed-ups  achievable  in  eigenvalue/eigenvector  computa- 
tions by  means  of  parallel  processing.  (Our  analysis  is 
applicable  to  Hermitian  matrices  since  they  can  be  con- 
verted to  real  symmetric  matrices  [11.)  There  are  two 
main  steps  involved  in  finding  the  eigenvalues/eigenvec- 
tors of  a symmetric  matrix: 

1)  Reduce  the  matrix  to  tri-diagonal  form  by 
Householder  transformations.  This  reqaires 

3 

approximately  2n  /3  multiplications  for  a 
nxn  matrix. 


2)  Apply  Francis-Kublanovskaja  QR  transforma- 
tions iteratively  to  diagonalize  the  tri- 
diagonal matrix  and  obtain  eigenvalues. 

Each  iteration  takes  3n  multiplications  and 
3n  divisions,  and  experiments  show  that  2 QR 
iterations  are  required  per  eigenvalue  [2]. 

First,  we  explain  Householder  transformations  briefly, 
and  show  how  the  computation  in  Step  1 can  be  speeded  up 
by  a factor  of  approximately  p,  the  number  of  processors. 

In  Step  1,  n-2  Householder  transformations  are  applied  and 
after  r-1  transformations,  we  have 
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Pr  is  selected  to  be  symmetric  and  orthogonal  and  in  the  prod- 

uct  PrAr-lPn’  the  circled  elements  in  A . are  annihilated. 
Note  if  r'y 


Householder  showed  that  for  any  given  column 
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There  exists  a symmetric,  orthogonal  matrix  P = 
such  that  w is  a nxl  vector  of  unit  length  (i.e 
and 


W.  W.  GflfFITMf-R 
RESEARCH  IMC. 


I-2ww 

, wTw  = 1 ) 


Pa 


* 

0 

0 


0 


l. 

2 2 2 2 
where  the  non-zero  element  is  (a,  + a^  + . . + a ) . 

'12  n ' 

Thus  Qr  can  be  selected  to  be  a Householder  transformation 
matrix.  (For  an  excellent  discussion  of  Householder  trans- 
formations, refer  to  Wilkinson  [1].) 


The  major  computation  involved  in  each  stage  is  the 
evaluation  of  the  product  QrBr_jQr-  Thus  we  consider: 

A = ( I-2wwT)B( I-2wwT) 


and  analyze  the  algorithmic  aspects  of  evaluating  A.  Let 
C = Bw  D = 2C-2w(CTw) 


Then  it  follows: 


A 


B- wDT 


If  A and  B are  mxm  matrices,  computation  of  C = Bw  and  wDT 

2 t h 

requires  2m  multiplications.  Hence,  in  the  r stage  there 

p 

are  approximately  2(n-r)  multiplications,  and  Step  1 requires 
3 

2n  /3  multiplications  on  a sequential  processor.  However,  in 
the  present  computer  with  p processors,  the  number  of  mtus 
required  to  perform  computation  in  the  rth  stage  is: 


3 
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and  hence  Step  1 requires: 
n-2 

2 prl  (n"r)  = 2n3/3P  + °(n2)  mtu 

r = l 

In  Step  2,  QR  transformations  are  applied  to  diagona- 

i.  L 

lize  the  tri-diagonal  matrix.  In  sL  stage  we  decompose 

As  - QSRS 

where  Q is  orthogonal  and  R is  upper-triangular,  and 
recombine  to  obtain 

As  + 1 ■ Vs 

since  R$  = Qs*As  = Q^A$,  we  have  As  + ^ = Q^ASQS.  This  pro- 
cess is  repeated  until  A$  converges  to  a diagonal  matrix. 
In  each  stage,  it  may  be  necessary  to  introduce  shifts  to 
accelerate  convergence,  as  in  the  following  fashion: 

W = Vs 

RsW  - As  + 1 

Ortega  and  Kaiser  [1]  gave  a method  for  simultaneous 
decomposition  and  recombination  in  QR  algorithm  for  tri- 
diagonal matrices.  Let 
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al  e2 
$2  «2 

. 0 

a 1 b2 
b2  a2  b3 
b3  a3 

0 

Then 

= “i 

pi  ■ n/ci-i  (if  ci-i  * °> 

= ci _2/0f  (if  cj-i  ” °) 

bi  ■ si-i<Pi  ♦ <*  * » 

si  " Biti/(Pi  * 8m):  ci  - Pt/(P,  ♦ «i+1) 

“i  " si f^i  + «i+1) 

‘i  = <i  * “i 

These  equations  are  executed  seriatim,  for  i + 1,2,...,  n, 
starting  with  cQ  = 1 and  uQ  = 0.  Also  the  value  of  3n  and 
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> 


n+1 


are  both  taken  as  zero. 


At  present,  there  is  no  known  method  for  exploiting 
parallelism  in  non-linear  recursion  equations  and  hence 

we  do  not  expect  a speed-up  of  p in  this  case.  However, 

3 2 

as  Step  1 requires  0(n  ) operations  as  compared  to  0(n  ) 

operations  in  Step  2,  the  overall  speed-up  appears  to  be 
primarily  decided  by  the  speed-up  of  Step  1.  Since  a 
speed-up  of  p is  achievable  in  the  first  step,  we  can  con- 
clude that,  in  general,  the  speed-up  for  obtaining  eigen- 
values is  p. 

A similar  analysis  applies  for  finding  eigenvectors. 
Since  matrix  multiplications  of  the  form 

(I-2wrwJ) (I-2wr  jwJj)  


are  involved,  we  expect  a speed-up  of  p in  the  first  step. 
In  the  second  step,  we  perform  multiplications: 

QT  0T  , • • • B 
ws  ”s - 1 


Thus,  using  the  techniques  presented  for  matrix  multiplica- 
tion, we  can  achieve  a speed-up  of  p. 
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[1]  Wilkinson,  J.  H.,  The  Algebraic  Eigenvalue  Problem, 
Clarendon  Press,  Oxford,  1965. 

[2]  Dahlquist,  G.,  et  al , Numerical  Methods,  Prentice-Hall, 
New  Jersey,  1974. 
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LISTING  OF  COMPUTER  PROGRAM  FOR 
MULTIPLE  SOURCE  LOCATION 


SI  GL4M.  F OR 


1 
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C MAIN  ROUTINE 

C LINK 

C SI  GL4M=SI  GLAM,  SYSLI  B/F/C 

C SI  GO!  1 /O:  1 /C 

C SI  GO  1 2/0: 1 /C 

C SI  GOl  3/0:  1/C 

C SI GOl 5/0: 1 

C SI  GOl  6/ SI  GOl  A/0:  1/C 

C SI  GO  1 7/0:  1 /C 

C SI  GO  18 /O:  1 

C 

IMPLICIT  DOUELE  PRECISION  ( A-H,  0-Z) 

LOGI  CAL*  1 TIMEIS?R) 

INTEGER  CASNUM/REINI  T,  UPDATE 
COMMON/ANSWER/  COMPUTC  10,  10)  , I C0MP1,  I COMPS 

COMMON/  SAVE  / DI  RVEC?  1 6,  1 6,  2)  / XM(  16/  16/2)/  SI  GVAR?  16)/  DEL  FAC 
1/  SI  GLOCC  1 to)  / WLIMI  T 

C0MM0N/C0MC0E/C0EFF?  16/  2)  / THRESH/ ANTNI  S 
COMMON/  JACO  /XMUC  32/  32)  / El  GEN?  32) 

COMMON/  CONSNT  / PI  / TWOPI  , ON EOFI 

COMM  ON  /H  2H  / CASN  UM  / I RND/ JRND/  KEINI  T,  UPDATE,  I TTY,  NUMANT/  I LOOP 
1/NUMSI  G/MSET, METHOD 
COMM  ON /NHL  TAP/ 1 HLNO,  IH  LN  2 
C 

CASNUM  = 0 
1 RND  = 2 69  49 
JRND  = -19A61 
PI  = 3.  14  1 59  2 6 5 358  9 79  D0 
TWOPI  = 2.  0D0*  PI 
ONEOPI  = 1 » 0D0/PI 
C TYPE  10 

C10  FORMAT?  ' INPUT  FILE  FOR  POLYNOMIAL  COMPUTATIONS  ON  THE  UNIT 

C 1 CIRCLE?’//) 

C CALL  ASSIGN? 3,  ’COEFF. FIL’/9) 

C 


M8E  is  bssi  qjattTT mcnoM 


TYPE  5000 

5000  FORMAT?’  TYPE  1 FOR  INPUT  FROM  TERMINAL  OR  2 FOR  INPUT  FROM 

1 FILE*  ’) 

ACCEPT~5001,  I TTY 

5001  FORMAT?  I 3) 

20  TYPE  5005 

50  0 5 FORMAT?////’  SELECT  ONE  OF  THE  FOLLOWING  ALGORITHMS:’,// 

1/’  TYPE  1 FOR  ALGORITHM  THAT  SELECTS  EIGENVECTORS 

2 Eased  on  eigenvalues.*// 

3/’  TYPE  2 FOR  ALGORITHM  THAT  SELECTS  F.I  GENVECTORS 

a Eased  on  inner  products.’,/ 

5, ’  TYPE  3 FOR  ALGORITHM  THAT  SELECTS  INDEPENDENT  VFXTORS  FROM’ 

6,  T MULTIPLE’, /,8X,  ’El  6ENVECT0R  SETS  BY  GRAM- S CM  I DT  PROCEDURE.' 

accept  5001 , Method  * 

c 

I FLAGS  1 
C 

I F?  I FLAG. NE.  1 ) GOTO  30 
CALL  SUBA 

30  I F?  I TTY.EQ.A)  GOTO  9990 


no  o n 


SIGL4M.F0R 
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I FC I TTY-EQ.2)  GOTO  40 
CALL  SUBBC  I FLAG) 

I FC  I FLAG.NE.  1)  GOTO  40 
CALL  SIGOPT 
GOTO  20 

40  CONTINUE 

CALL  SI  GO  1 2 


I FC  METHOD.  EQ.  1)  CALL  SIG018 
I K METHOD.  EQ.  2)  CALL  SIG017 
I FC  METHOD.  EQ.3)  CALL  SI  GO  I 6 


CALL  SUBCI  KC  NUMANT-  1 ) 

GOTO  20 

9990  CONTINUE 
STOP 
EM  D 

BLOCK  DATA 

COMM OM/MRLTAP/I  FILNO,  I FILN2 
DATA  I FILMO/- 1/ 

END 

SUBROUTINE  CMPLXMC  A,  B,  C,  D,  K,  F ) 

DOUBLE  PRLCI  SI  ON  A,  B,  C,  D,  E,  F 

E=  A*  C-  B*  D 

F=  A*  D+  B*  C 

RETURN 

END 
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SUBROUTINE  SURA 

IMPLICIT  LOUELE  FI  EC  I SI  ONC  A-H,  0-£  ) 

REAL  ElOAT 

INTEGER  REI  MI  T,  UPDATE,  CASNUM 

COMMON /ANSWER/  COMPUTC  1 0,  1 0)  , I C0MP1  , I COM Pg 

COMMON/  SAUL  / DI  PVECC  1 6,  1 6,  2)  ,XMC  1 fc,  1 6/  2)  , SI  GVARC  1 6)  , DELE  AC 
1,SI  GLOCC  16),  VLIMIT 

COM  MOM/  COMCOE  / COLE  K 1 6,  2)  , THRESH,  ANT-MI  S 
COMMGM/JACO/XMUC  32,  32)  , El  GE'MC  32) 

COMMON!/  COM SM T / FI  , TVOPI  , ONEOPI 
COMMON/HH/  AIMC  1 6,  2)  , SI  G(  1 6,  2)  , SMC  1 6,  2) 

1#  XI  C 16,  16,  2) 

COMMON! /H24/CASMUM,  I 1ND,  JHND,  REI  MI  T,  UPDAl  E.  I TTY  , MUM  AMT,  I LOOP 
1,  MUM  SI  G,  M SET,  METHOD 
COMMON /iMRLTAP/I  EILMO,  I EILM2 
C 

GOTO  (1000,1200,1300,9990)  ITIY 


1000  COMTIMUE 

TYPE  3000 

30  0 0 EORMATC /////,  IX,  ’NUMBER  OE  AMTEMMA  ELEMENTS?') 

ACCEPT  300 1 , MUMAMT 

3001  EORMAT< 17) 

C 

C SET  EOR  MULTIPLE  ANALYSIS 

C 

I E(  METHOD.  GT.  1 )NUMAMT=NUMAMT+  1 
I E(  MUMAMT.  GT.  IS)  GOTO  1000 
C 

TYPE  3002 

3002  E0RMAT(  IX,  ’DELAY  EACTOR  2*  D/LAM  BDA  ?’) 

ACCEPT  3010, DELE  AC 

3010  EORMATC  E 1 A.  7) 

TYPE  3003 

3003  EORMATC  IX,  ’MUMPER  OE  SIGNAL  SOURCES?’) 

ACCEPT  3001 , MUM SI G v 

TYPE  3004 

3004  EORMATC  IX,  ’SI  GMAL  SOURCE  LOCATI  OMS  THETA/PI?’) 

ACCEPT  3005, ( SI  GLOCC  I ) , 1 = 1,NUMSI  G) 

30  05  EORMATC  5C  E 14.  7)  ) 

TYPE  3006 

3006  EORMATC  IX,  ’SI  GMAL  SOURCE  V/ARI  AM CES ? ’ ) 

ACCEPT  3005,  C SI  GVARC  I ),  I = 1,NUMSI  G) 

TYPE  3007 

3007  EORMATC  IX,  ’UPDATE  INTERVAL  IN  NUMBER  OE  SAMPLES?’) 

ACCEPT  3001,  UPDATE 

TYPE  3003 

3003  EORMATC  IX,  'REINITIALIZATION  INTERVAL  IN  NUMBER  OE  SAMPLES?’) 
ACCEPT  3001,  REI  NIT 
TYPE  3307 

3307  EORMATC*  INTERNAL  AMTEMMA  NOISE  EACTOR?’) 

ACCEPT  3010,  AMTNIS 
I THRSS=  2 

I EC  METHOD.  EQ.  2)  GOTO  3530 


SIG01  1 • E OR 


4 


8 6- OCT- 7 7 13:  07:21  PAGE,! 

( 

THIS  PACE  IS  BEST  QUALITY  PRACTICABLE 
FftQB  CWPY  FUR3U.SHED  TO  DOC ' 

3505  TYPE  50550 

ACCEPT  50551, ITHRSS 

3630  I H IT4RSS.LT.  1.  OR.  ITHRSS. GT. 2)  GOTO  3505 
IK  I THESS.  EC.  1)  THRESH=0.  0D0 
IK  ITHRSS.  EQ.2)  TYPE  50510 

C3308  EORMATC  * THRESHOLD  - ANTENNA  INTERNAL  NOI  SE  POWER?’) 

I EC I THRSS. EQ. 2) ACCEPT  3010, THRESH 
I HMETHOD.LT.  3)  GOTO  3530 
TYPE  3330 

3330  EORMATC /, * TOL ERANCE  E OR  GRAM-SCHMI DT  ORTHOGONAL  I ZATI ON ’, 
l*  PROCEDURE  ?’> 

ACCEPT  3010, WLIMIT 

3530  CASNUM  = CASNUM+l 
DELE  AC  = DELEAC*  PI 
GOTO  1 400 


Cl  200  EX  PAM  SI  OM  EOR  El  LE  INPUT. 

1200  CONTINUE 

1300  CONTINUE 
C 

MUM  AN  T=  A 
DELEAC= 1 • 0D0 
DELEAC=  deleac+pi 

I E(  I EILMO.  GT-  0)  GOTO  50500 
13010  TYPE  1301 

1301  EORMATC  ////,  ’ MAG  TAPE  EILE  NUMBERS  EP.OM  ...  TO  ...?’) 

ACCEPT  50001,1  EILN  1,1  El LiM 2 

50001  EORMATC  21  5) 

50549  TYPE  50550 

50550  EORMATC’  TYPE  l EOR  INTERACTIVE  THRESHOLD  SELECTION.’,/ 
1,’  TYPE  2 EOR  NON-INTERACTIVE  THRESHOLD  SELECTION.’) 
ACCEPT  50551, I THRSS 

505  51  EORMATC  ID 

I EC  I THRSS. LE.  0.  OR.  I THRSS.  GT.  2)  GOTO  50549 
I EC  I THRSS.  EQ.  1)  THRESH=0.0D0 
I EC  I THRSS.  EQ.  2)  TYPE  50510 

5051  0 EORMATC’  THRESHOLD  EOR  SELECTING  EIGENVECTORS?’) 

I EC  I THRSS.  EQ.2)  ACCEPT  5051  1,  THRESH 
5051  1 EORMATC  D14.  7) 

I EILN 0=1 EILN 1-1 
50500  I EILNO=I EILNO+1 

I EC  I EILNO.  GT.I  EILN2)  GOTO  13010 
TYPE  50003, I EILNO 

50003  EORMATC///,’  NRL  MAG  TAPE  EILE  NUMBER  : * , I 2) 

C 

TYPE  50004, NUMANT,  DELEAC/PI 

50004  EORMATC///,’  ANTENNA  ELEMENTS  : ’ ,11,/, 

1 * DEJLFAC  2*  D/LAMBDA  :’  ,1PDA.1,/, 

2 r NUMBER  OE  SAMPLES  : 1024’) 

CALL  SIGMTR 

RETURN 

C 

1400  CONTINUE 

I EC  REINI  T/UPLATE  .LE.  10)  GOTO  1420 


SI  GDI  1 « EQR 
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1420 
1440 
40  0 2 

4003 

40  04 
4P0b 

4006 

C 

1500 

4 

5 

C 

C 


10 


999  0 


STOP  'TOO  MANY  STEPS' 

I ECNUMSI  G-LE.  10)  GOTO  1440 
STOP  'TOO  MAM  If  SI  UMAL  SOURCES ' 
CONTINUE 

WRI  TEC  2,  4002)  CASNUM 
ECRMATC  1H1, ///,  ' SIMULATIOM  OE 
1 PROBLEM  ', ///,  ' CASK  MUMPER  : 
WRI  TE(  2,  4003)  MUM  AMT 
EORMATC ' AMTEMMA  PARAMETERS  S ' 
112) 

WRI  TEC  2,  4004)  DEL  E AC* ON EOPI 
EORMATC'  DELAY  EACTOR  =’ 

WPI  TEC  2.  4005) 

EORMATC  ' THE  COVARIANCE  MATRIX 
1 NOISES  IS  ASiUMED  TO  BE',/,' 
WRI  TEC  2.  4006)  , THRESH 
EORMATC///,'  ALGORITHM  PARAMETER 
1 AMTEMMA  I MTERMAL  MOISE  POWER  = 


TS  BBSt  QiiKtt  n mcncASu 


MULTI  - SOURCE. 

',  I 5/  ///) 


DIRECTION  EIMDIMG 


//, 


NUMBER  OE  ELEMENTS 


= ’,  1PD10.3) 


OE 

AM 


I MTERMAL 
I DENT  I TY 


AMTE-MMA 
MATRIX • 


, //,  3X,  'THRESHOLD 
1 PD10.  3) 


I LOOP  = 0 
I COMPl=  0 
ICOMPP=0 

DO  4 1 = 1#  NUMSIG 
SIGLOCCI)  = PI  *SI  GLOCC  I ) 

DO  5 I = 1,  MUMAMT 
M b J = 1*  MUMAMT 
DO  5 K = 1,2 
XMC  I , J,  K)  = O.0D0 

MODULE.  DIR-VEC-GEM.  GEMERATES  DIRECTION  VECTORS  OE  SIGNAL 
SOURCES. 

DO  10  I DIR  a J,  NUMSIG 

\X  = -DELEAC*DSIMC  SI  GLOCC  I DI  R)  ) 

DO  10  JL'IR  = 1,  MUMAMT 
XY  = DEE  EC  ELOATC  JDI  R-  1 ) )*XX 
DIRVECC  I DI  1,  JDIR,  1)  = DCOSCXY) 

DI  RVECC  I DI  I,  JDI  R,  2)  = DSIMCXY) 

CONTINUE 

MELAG=»  0 
RETURN 
I TTY=»4 
RETURN 
END 

SUBROUTINE  ‘UFBC  I ELAG) 

IMPLICIT  DOUBLE  PRECI  SI  ONC  A-H,  O-t;  ) 

REAL  RAN, ELOAT 

INTEGER  REI N I T,  UPDATE,  CASNUM 

COMMON /COMCOE/  COEEEC  1 6,  2)  , THRESH,  AN  TNI  S 

COMMON/  SAVE  / DI  RVECC  1 6,  1 6,  2)  , XMC  1 6,  l 6,  2)  , SI  GVARC  1 6)  , DEL  E AC 
1 , SI  GLOCC  16) 

COMMON/ JACO/XMUC  32,  32)  , El  GENC  32) 

COMMON/  CONSNT  / PI  , TWOPI  , ON EOpI 
COMMON/HH/  AINC  16,  2)»SIGC  16,  2),SNC  16,  2) 

1*  XI  Cl  6,  16,  2) 


SI  601  1.F0R 
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C0MM0N/H2H/CASN  UM,  I RND,  JRND,  RLHIT,  UPDATE,  I TTY,  NUMANT,  ILOOP' 

1 , N UM  S I G 
C 

I FLAG=0 

20  00  ILOOP  = I LOOP  + 1 

C NODULE  ANT- 1 NT-NOI SE.  GENERATES  ANTENNA  INTERNAL  NOISE  SAMFLES. 

DO  20  I ANT  = 1,  NUMANT 
ANTI  = DELECRANC  I END,  JRND)  > 

AN  T 2 = DEL  E(  R AN  C I RN  D,  J RN  D)  ) 

CALL  GAUS(  ANTI,  ANT2,  AIN(  I ANT,  1 ) , AI  N(  I ANT,  2)  , ANTNI  S) 

20  CONTINUE 

C MODULE  SI  G- GEN.  GENERATES  RECEIVED  SIGNAL  SAMPLES. 

DO  30  I SI  G = 1,  NUMSI  G 

SI  G 1 = DEL  EC  RANI  I RN  D,  J RN  D)  ) 

SI  G2  = DEL  EC  RANC  I RN  D,  JKN  D)  ) 

CALL  GAUSC  SI  Gl,  SI  G2,  SI  G(  I SI  G,  1)  , SI  6<  I SI  G,  2)  , 0.  bD0* 

1SI  GVARC  I SI  G)  ) 

30  CONTINUE 

DO  40  I SI  G = 1,  NUMANT 
SNCISIG,  1)  = 0.0D0 
SNCISIG,  2>  = 0.0D0 
DO  40  JSI  G = 1,  NUMSI  G 

CALL  CMPLXMC  SI  GC  JSI  G,  1 ) , SI  GC  JSI  G,  2)  , DI  RVECC  JSI  G,  I SI  G,  1)  , 

1DI  IVECC  JSIG,  ISIG»2),A1,A2> 

CALL  CMFLXAC  SNC  I SI  G,  1 ),  SNC  I SI  G,  2),  Al,  A2,  SMC  I SI  G,  1 ) , 

1 SNCISIG,  2)) 

40  CONTINUE 

C MODULE  COMBINE.  COMBINES  ANTENNA  INTENAL  NOISE  AND  RECEIVED 

C SIGNALS. 

DO  50  I COM  = 1,  NUMANT 

SNC  I COM,  1)  = AINCI  COM,  1>  ♦ SNC  I COM,  1) 

SNC I COM, 2)  = AINCI COM, 2)  + SNC  I COM, 2) 

50  CONTINUE 

C MODULE  COM-COVAR-MAT.  COMPUTES  THE  COVARIANCE  MATRIX  01-  THE 

C COMBINED  SIGNAL  SAMPLES.  ( 

DO  60  I COM  = 1,  NUMANT 
DO  60  JCOM  = 1,  NUMANT 

CALL  CMPLXMC  SNC I COM,  1 ), -SNC I COM, 2> , SNC  JCOM,  1 ), SNC  JCOM, 2) , 

1A1, A2) 

XMC I COM, JCOM, 1 ) = XMC I COM, JCOM,  1 ) + Al 
XMC  I COM,  JCOM,  2)  = XMC  I COM,  JCOM,  2)  + A2 

60  CONTINUE 

C MODULE  UPDATE-COVAR-MAT.  UPDATES  COVARIANCE  MATRIX  OH 

C REINITIALIZES  TO  ZERO  MATRIX. 

I*  CMODC  ILOOP,  FEINI  T)  .NE.  0)  GOTO  3000 

I PLAG= 1 

RETURN 

3000  CONTINUE 

IF  CMODC  ILOOP,  UPDATE)  .NE.  0)  GOTO  2000 
DO  70  I UPD  = 1,  2* NUMANT 


SI  liOl  1.  RUR 
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70 

D 

7 b 


80 


XHIS  PiOE  IS  BEST  «U*LItI  mCHUS* 

rooa  oopt  jvjRtkisHBu  ro  doc  — *■" 

DO  70  JUPD  = 1,  2*N  UMAN T 
XMUC  I UPD,  JUPD)  = 0.0D0 
DRLOOp=  DBLEC  RLOATC  ILOOP)  ) 

TYPE  7 b,  ( < ( X M ( I UPD,  JUF'D,  KUPD)  / DRLOOP,  KUPD=  1 , 2;  , JUPD.=  1 , NUMANT)  , 
II  UPD=  1,  NUMANT) 

R ORM  ATC  8 ( IX,  1-7.  A)  ) 

DO  80  IUPD  = 1,  NUMANT 
DO  80  JUPD  = 1,  MUM  AMT 

XMUC  I UPD,  JUPD)  = XMC  I UPD,  JUPD,  1 ) /DRLOOP 

XMUC  I UPD+NUMANT,  JUPD+NUMANT)  =XMC  I UPD,  JUPD,  1 )/ DRLOOP 

XMUC  I UPD,  JUPD+NUMANT)  = -XMC  I UPD,  JUPD,  2)  / DI-LOOP 

CONTINUE 

I RLAG=0 

RETURN 

END 

SUBPOUTINE  CMPLXAC  A,  E,  C,  D,  E,  R > 

DOUBLE  PRECISION  A,  B,  C,  D,  E,  R 

E=  A+  C 

R = B+D 

RETURN 

END 

SUBF.OUTI  N E OMPLX  DC  A,  B,  C,  D,  E,  R ) 

DOUBLE  PRECISION  A,  B,  C,  D,  E,  R,  G 

G=  C*  C+  D*  D 

E=A+C+B*  D 

R = B*  C- A*  D 

E=E/G 

R=  R /G 

RETURN 


END 

SUBROUTINE  CMPLXTC  A,  B,  C,  D) 

DOUBLE  PRECISION  A,  B,  C,  D 
C=  DSQRTC  A*A+F*P) 

D=  DAT  AN  2C  B,  A) 

RETURN 

END 

SUBROUTINE  GAUSC  Al,  A2,  Bl,  82,  SI  G) 
IMPLICIT  DOUBLE  PRECI  SI  ONC  A-B,  O-Z  ) 
COMMON/  CONSNT  / PI  , TV OPI  , ON EOPI 
A=  - 2*  0D0*  SI  G*  DLOGC  Al) 

A= DSQRTC  A) 

Bl  = A*  DC 0 SC  TU0PI  + A2) 

B2  = A*  DSINC  TV.0PI*A2) 

RETURN 

END 

SUBROUTINE  CMFLX  CC  A,  B,  C,  D) 

DOUBLE  PRECISION  A,  P,  C,  D 
C=  A*  DCOSC  B) 

D=  A*  DSI  NC  B) 

RETURN 

END 


THIS  PACE  IS  BEST  QUALITY  PRACIICABil 
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SUBROUTINE  SI  GO 12 

IMPLICIT  DOUBLE  PRECISION  CA-M,0-t:> 

REAL  FLOAT 
C 

^ COMMON /SAVE/ DIR  VECC  1 6,  1 6,  2>,XMC  1 6,  1 6*  2>,  SI  GVAKC  1 6)  , DEL  FAC 

l»  SIGLOCC  16) 

COMMON/  JACO/XMUC  32*  32)  » El  GENC  32) 

C 

INTEGER  CASNUM,REINI  T,  UPDATE 

C0MM0N/H2H/CASNUM,  I RND,  JHND,  REI  NIT,  UPDATE,  I TTY  , NUMAMT 
1,  ILOOP,NUMSl  G,MSET, METHOD 
C 

C UFE  XMU2  TO  GEN  El  GLNP  VIA  JACOBI 

C 

CALL  JACOBI  < 2*NUMANT,  50,  1 . 0D-  10,  1 . PD-  IP,  1 • 0D-  1 0) 

I FC  METHOD.  EQ.  1)  RETURN 
CALL  ASSI GN(  1 , ' SI GP2»  TMP ' , 9) 

WRI  TE<  1,  1 0)NUMANT 
10  F0RMATCI3) 

WRI  TEC  1,  l 5)  CC  XMUC  I , J)  , J=  1 , 2*NUMANT>  ,1  = 1, 2+NUMAMT) 

15  FORMATC  IX,  ADI  5.8) 

WRI  TEC  1,  1 5)C  El  GENC  I ),  1=  1 , 2*NUMANT) 

CALL  CL  OS  EC  1 ) 

C 

C GEN  TASS  l DATA  NOW 

>C 

DFL=  DBLEC  FLOATC  I LOOP)  ) 

NMAN  T 1 =.N UM AN T-  J 
DO  8 0 1 = 1,  NMANT  1 
DO  8 0 J=  1 , NMAN  T 1 
XMUC  1 , J)  =XMC  I , J,  1 ) /DFL 
XMUC  I+NMANTl,  J+NMANT 1 ) =XMC  I , J,  1 ) / DFL 
XMUC  I , J+NMANTl  )=-XMC  I , J,  2)  / DFL 
80  CONTINUE 

C 

CALL  JACOBI  C 2+NMANT1, 50,  1 . 0D-  10,  1 . 0D-  IP,  1. 0D-  10) 

RETURN 

END 

SUBROUTINE  JACOBI CN, I TM AX, EPS  1 , EPS2,  EPS3) 

IMPLICIT  DOUBLE  PRECI  SI  ONC  A-H  , 0- Z ) 

COMMON /JACO/AC  32,  32)  , El  GENC  32) 

DIMENSION  TC  32,  32),AIKC  32)>N0C  32) 

LOGICAL  FINE 
C TYPE  6666,  A 

D6666  FORMATC  * JACO  ' , /C  5C  1 PE9 . 3)  ) ) 

DO  A I = 1 , N ‘ 

DO  A J = 1,  N 
A TC  I , J)  = 0.0D0 

NMI=N-1 
SI  GMA 1 = 0»  0 
OF  F DSQ=  0*  0 
DO  5 I = 1 , N 

SI  GMA  1 = SI  GMA l ♦ AC  I , I )**2 

TC  I , I ) = 1 . 0 

IP1-I+I 


V 


?' — ’r-'T_ 


SI G012. 1-OR 
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DO  5 J=  I P 1 * N yxM  (x»PY  rURAISHED  TO  DUC  — - 

b OH-DSQ=OH-DSQ+AC  I*  J>**2 

6 S=2«  0*  01*  h DSQ+SI  GMA1 

DO  26  I TER= 1*1 TMAX 
DO  20  1 = 1*  MM  1 
IP1  = I + 1 
DO  20  J= I P 1 * N 
Q=  DABSC  A< I*I)-AC J* J>> 

I > < Q-LE.EPS1)  GOTO  9 

II-  C DABSC  AC  I * J)  ) .LE.  EPS2)  GOTO  20 

P=2.  0*AC  I,  J)*0/C  AC  I*  I )-AC  J*  J)  ) 

SPQ=  DSQRTC P* P+G*  Q) 

CSA= DSQRTC  ( 1 . 0+G/SPQ) /2«  0) 

SMA=P/C  2. 0*  CSA*  SFQ  ) 

GOTO  10 

9 CSA= 1 . 0/DSQRTC 2. 0DO) 

SNA=CSA 

10  CONTINUE 

DO  11  K= 1 * N 
HOLDKI  = TC  K>  I ) 

TCK*  I ) =HOL  DK I * CSA+  TC  K*  J)  * SNA 

11  TC  K*  J)  =HOLDKI  * SNA-  TC  K*  J) * CSA 
DO  16  K=I*N 

It  CK.GT.J)  GOTO  IS 
AI  KC  K)  = AC  I * K ) 

AC  I * K)  = CSA*  A I KC  K ) + SN  A*  AC  K*  J) 

If  CK.NE.J)  GOTO  1 A 

AC  J*K)  = SNA*AIKCK)-CSA*AC  J*K> 

14  GOTO  16 

15  HOLDI  K=  AC  I * K) 

AC  I * K)  = CSA*HOLDI  K+  SNA*  AC  J*  K> 

AC  J*  K>  = SNA*HOLDI  K- CSA*  AC  J*  K) 

16  CONTINUE 

AI  KC  J ) = SNA* AI  KC  I > - CSA*  AI  KC  J) 

DO  19  K=  1 * J 

If  CK.LE.I)  GOTO  18 

AC  K*  J)  = SN  A*  AI  KC  K)  - CSA*  AC  K*  J) 

GOTO  19 

18  HOLDKI  = AC  K*  I) 

ACK*  I ) = CSA* HOLDKI  + SNA*  AC  K*  J) 

AC  K*  J)  = SN  A*H  OLDKI  - CSA*  AC  K*  J) 

19  CONTINUE 

80  ACI*J>  = 0.  0 

SIGMA2=0.  0 
DO  21  I = 1 * N 
EIGENC  I )=AC  I,  I > 

21  SIGMA2=SI  GMA2+EI  GENC  I >**2 

IP  C 1.0-SIGMA1/SIGNA2.LT.EPS3)  GOTO  100 
26  SI  GMA1  = SI  GMA2 

TYPE  200 

200  fORMATC  IX*  'CONVERGENCE  DID  NOT  OCCUR’) 

RETURN 

C PROGRAM  TO  SORT  THE  EIGENVALUES  IN  DESCENDING  ORDER. 

100  DO  500  I = 1*  N 


SI  GO  1 2.  F OR 
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N OC I > = I 
FINE  = . TRUE. 

00  bl0  J = 2,  N 

IF  < EIGENC  J-  1>.GE.  EIGENC  J)  ) 

FINE  = . FALSE. 

TMP  = EIGEN(J-l) 

EIGEMCJ-1)  = El  GENC  J) 

EIGENC  J)  = TMP 
NTMP  = NOC  J-  1) 

NOCJ-1)  = NOC  J) 

NOCJ)  = NTMP 
CONTINUE 

IF  C .NOT.  FINE)  GOTO  501 
TYPE  700, N/2,  C El  GENC  I >,  I = 1 
FORMAT?  //,  ' EIGENVALUES  FOR 
1,C  IX, 5 F 1 4 • 7 ) ) 

DO  650  I = 1,  N 
DO  650  J = 1,  N 
AC  J,  I ) = TCJ,NO(  I ) ) 

RETURN 

END 


1HIS  PACE  IS  BEST  QUALITY  PRACTICABLE 
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GOTO  510 


ANTENNA  ELEMENTS 


SIG013.F0R 


C 

D 

E5500 


10 


200 

210 
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ISIS  PA&E 
FROM  OOP* 


B BIST  W1U.IH 
yURUlSHED  10  DDC  "" 


SUBROUTINE  SUBCIR(N) 

IMPLICIT  DOUELE  PRECISION  !A-H,0-Z) 

REAL  FI„  OAT 

INTEGER  CASNUM,REI  NIT#  UPDATE 
COMMON/ANSWER/  COMPUT!  1 0,  1 0)  , I COMP  1 » I C0MP2 
COMMON/  COMCOE  / COEF  F!  1 6,  2)  , THRESH 
COMMON/  CONSNT  / PI  , TVOFI  , ONEOPI 

COMMON/SAUE/  DI  RVEC!  1 6,  1 6,  2),XM!  1 6,  1 6,  2)  , SI  GUAR!  1 6),  DEL  FAC 
1,SIGL0C!  16),  VLI  MIT 

C OMM  ON  /H  2H  / C ASM  UM  , I RMD,  JRMD,  REINI  T,  UFDATE,  I TTY,  NUMANT,  I LOOP 
1,NUMSI  G,MSET 
REAL  A(  1024,2) 

DIMENSION  REILT!  15),  PHI  VAL!  1 5) 


TYPE  5500, MSET 

FORMAT!  * SUECIR,  MSET=',I3) 

IN  DX  = 0 * 

I COMP 1 = 1 C0MP1+ 1 
I C0MP2=  0 

TPI  DNU=  TVJOPI  / 1 024 • 

CALL  ASSIGN!  3,  * COEFF. FI L’, 9) 
REWIND  3 
READ!  3)  A 
CALL’  CLOSE!  3) 

X 1 = DEL  E!  A!  1024,  1)  ) 

X2=DBLE!  A!  1084,2)) 

RESSAV=  CAL PUN! X 1 ,X  2, N) 

SLPSAV=  1 • 0D0 

DO  300  LOOP  = 1,  1024 

XI  = DELE  ! A!  LOOP,  D) 

X 2 = DELE  ! A!  LOOP,  2)) 

R1  = COEFHN+1,1) 

R2  = COEF  F(  N+  1,2) 

DO  10  I = N,  1,  -1 
Rll  = R1*X  1 - R2*X2 
R12  = R1*X2  + R2*X 1 
R1  = COEFF! I, 1)  + Rll 
R2  = COEFF!  1,2)  + R12 
CONTINUE 

RE  i=  R1*R1  + R2*R2 
SLOPE=RES-RESSAV 
IFILOOP.EQ.  1)  SLPFST=  SLOPE 
I F!  SLOPE.LT.0.0D0)  GOTO  290 
I F! SLOPE. EQ. 0. 0D0)  GOTO  200 
I F! SLPSAV. GE. 0. 0D0)  GOTO  290 
DFL=DBLE!  FLOAT!  L OOP- 1 ) > 

DFLM2=  DFL-2*  0D0 
CONTINUE 
X 1=  TPI  DNU*  DFLM2 
X2=  TPI  DMU+DFL 
INDXMNIK+1 

CALL  FE0NAC!X1,X2,  RESLT!  INDX)) 
X 1=  DCOS!  RESLT!  INDX)) 

X2=  DSI  N!  RESLT!  IN  DC)  ) 

PH  I A=  C AL  F UN  ! X 1 , X 2,  N ) 


SI  G013.  FOR 


D 

L6066 

C 


240 

C 

C 

245 

251 

D 

D6029 

90 

300 


31  0 
C 

339  0 


C 

C 


345 


340 

C 
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TYPE  6066#  PHI  A>  RESLTC  I N DX ) 

FORMAT?  ' PHI  A=  ’#  1PD16.8#  * RESLT=  ’ # D1  6*  8 ) 

FHI  VAL?  I N DX  > = PH  I A 
I F?  INEK.LE.MSET)  GOTO  290 
DO  240  1 = 1 # MSET 

I F? PHI A.GE.PHI VAL?  I > > GOTO  240 

PHI  VALC  I ) = PHI  A ! FOUND  A LESSER  VALUE 

RESLTC  I > = RESLT?  INDX) 

GOTO  245 
CONTINUE 


INDX  = INDX-  1 
GOTO  29  0 
CONTINUE 

TYPE  6029#  RESLTC  INDX5#  PHI  A 

FORMAT?  ’ A= ’#  1PD1  6-9#  ’#  PHI  OF  A='#D15.9) 

SLPSAV=  SLOPE 

RESSAV=HES 

CONTINUE 

1 F?  SLOPE.  GE.  p.  0)  GOTO  310 

I F?  SLPFST.LE.  0.  0)  GOTO  310 

DFL=0.  0 

DFLM2=-2.  0D0 

GOTO  210 

CONTINUE 

TYPE  339  0 

FORMAT?//#’  ESTIMATES  OF  ANGULAR  LOCATIONS  :’) 

DO  340  I = 1 #T  N DX 

Al  = RESLTC  I ) / DEL FAC 

IF?  Al.GT.  1.OD0)  A 1 = A 1 - 2.  0D0 

A 1 = DAT  AN  ? Al/DSQRT?  1 . 0D0- A 1 * A 1 ) > * ONEOPI 

I C0MP2=I C0MP2+ 1 

CQMPUT?  I C0MF2#  I C0MP1  )=A1 

SQRPHI  = DSQRT? PHI VAL? I ) ) 

TYPE  345#I»A1»  SQRPHI 

FORMAT? /#  5X#  'SOURCE  ’,12#’  : ’#  / 

\»  5X#  * ESTIMATED  ANGULAR  LOCATION  IN  RADIANS/PI  :’#F11.7#/ 

2#  5)?#’  RESIDUAL  ERROR  ON  UNIT  CIRCLE  :',1PD11.4> 

CONTINUE  * 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  CAL  FUN?  X 1 #X2#  N) 

IMPLICIT  DOUBLE  PRECISION  CA-H.O-Z) 

COMMON/  COMCOE  / COEF  F ? 1 6#  2)  # THRESH 
Rl  = COEFFCN+1#  1) 

R2  = COEF F?  N+  1 # 2) 

DO  10  I = N#  1#  -1 
Rl 1 = R1*X1  - R2 *X2 
Rl 2 = R1*X2  «■  R2*X 1 
Rl  = COEF  F? I # 1)  ♦ Rl 1 
K2  = COEF  F?  I,  2)  + R12 


i8no§  o o ono  n o n o c o o 


SI  G013. FOR 


I 
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CONTINUE 

CAL  F UN  = Rl  + Rl  + R2*R2 

RETURN 

END 


IRIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FRO*  OOPY  YUKWSHEC  TO  DSC 


SUBROUTINE  F BON  ACC  AIN,  BI  N,  RESLT) 


FIBONACCI  SECUEN  CE 


INFUTS  AIN,  BIN 
OUTPUT  RESLT 


IMPLICIT  DOUBLE  PRECISION  CA-H,0-Z) 

INTEGER  CASNUM,  REIN  I T,  UPDATE 
COMMON/  COMCOE  / COEFFC  16,  2),  THRESH 

C0MM0N/H2H/CASNUM,  I RND,  JRN  D,  REI  NI  T,  UPDATE,  I TTY,  NUM ANT,  I LOOP 
1 , N UM  S I G 

DIMENSION  T 1 ( 60) , T2C  60)  , TAUC  60)  , AC  60),  E<  60) 

DATA  IRSTFL/0/ 

DATA  MAX  / 60/ 

COMPUTE  FIBONACCI  SEQUENCE 

I PC  I RSTFL.NE.0)  GOTO  29 
I RSTPL=  1 

TAUC 1 ) = 1 • 0D0 
TAUC  2)  = 1 • 0D0 

DO  20  K=  3,  MAX  ! COMPUTE  TAU 

TAUC K)= TAUC K-  D + TAUCK-2) 


AC  1 ) = AI  N 
BC  i)  = BIN 
C 

N=N  UM  AN  T-  1 
C 

DO  40  K=l, MAX-2 

T1CK+  1)=C  TAUC  MAX-  1-K)  /TAUC  MAX*  1 -K)  )*  C BC  K)  - AC  K)  ) +AC  K) 
T2CK+  1)=C  TAUC  MAX-K)  /TAUC M AX+ 1 -K) ) *C  BC  K)  - AC  K)  ) + AC  K) 

C 

X 1 = DCOSC  TIC  K+ 1 ) ) 

X2=DSINC  TlCK+1)) 

R 1 = COEFFCN+1,1) 

R2  = COEFFCN+1,2) 

DO  10  I = M,  1,  -1 
Rll  = Rl*Xl  - R2*X2 
R1 2 = R1+X2  ♦ R2*X 1 
R1  = COEFFC I,  1)  + Rll 
R2  « COEFFC 1,2)  ♦ R12 
10  CONTINUE 

PHI  1=  R1*R1  + R2*R2 
C 

Xl=  DCOSC  T2C  K+- 1 ) ) 


CM 


I 


SI  G014.  hOh 


26-0CT-/7  13:48:14  PAGE 
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C 

C 

c 


c 

c 

c 

D 

Lb 

D 

D6 


617  0 


618  0 
D 

D618  1 

618  5 

C 

C 

C 


619  0 


6199 

C 

D 


SUBROUTINE  SI  GABT 


rua  is  best  in  moucUi* 


ABOVE  THRESHOLD  El  OEM  SELECT  I OM 


IMPLICIT  DOUBLE  PRECI  SI  0M<  A*H»  0-7,) 

REAL  FLOAT 

INTEGER  RE1NI  T,  UPDATE,  CASNUM 
COMMON/ COM  COE/  COEhK  16,  2),  THRESH 
COMMON/  CO\JSNT  / PI  , TWOPI  , ONEOPI 
COMMON /0 16/  AINC  16,  2),SIGC  16,2),SNC  16,2) 

1,  XI  C 16,  16,  2) 

COMMON /H2H/CASNUM,  I HMD,  JRND,  REINI  T,  UPDATE,  I T’lY,  NUMANT,  I LOOP 
i,NUMSI  G, MSET 

COMPUTE  LAMBDA,  1UM  Oh  THE  SGUARES  Oh  W 


TYPE  S,  NUMANT,  MSET 

hORMAT< ' SI  GABT-NUMANT-  * , I 3,  ’ MSET=’,I3) 

TYPE  6,  C ( < XI  < I > J,K),K= 1, 2) , J= 1,MSET)  , I = 1,  NUMANT) 

FORM  AT<  • SI  GABT-XI  =*,/,C  1X,4D1S.R>) 

DO  6 1 70  I = 1 , N UM AN  T-  1 
DO  6170  J= 1,MSET 
XI< I , J, 2)=-XI< I , J, 2) 

CONTINUE 
SUMSQV=0.  0D0 
I = NUMANT 

DO  6180  J-  1 , MSET 

SUMSQW-SUMSGW+XI ( I , J, 1 > *XI C I , J, 1 ) +XI ( I , J, 2) *X I C I , J, 2) 

CONTINUE 

TYPE  618  1 , SUMSQW 

hORMATC  ’ SUMSQW= ’,  1PD20.  12) 

I H DABSfl.ODO- SUMSQW). LT.1.0D-S)  TYPE  618  S 

hORMAT</,  ’ WARNING  ABSOLUTE  VALUE  Oh  l-SUMSGW  IS  LESS  THAN  1.0D- 
1S‘> 

COMPUTE  COEFh 

COEh  K NUMANT,  1)=  1.0D0 
Ch«=-( 1.0D0- SUMSQW) 

COEh  K NUMANT,  2)  «=0.  0D0 

MP1-NUMANT 

DO  6199  I o 1,  NUMANT-  1 

COEh  K I , 1 ) =0.  0D0 

COEhK  I , 2)  = 0«  0D0 

DO  6190  1,MSET 

COEh  h<  I , 1 ) >=  COEhK  I , 1 ) +XI  ( MP1 , J,  1 ) *X  I < I , J,  1 ) 

1 -XKMP1,  J,  2)*XI  ( I ,<J,  2) 

COEhK  1,2)  = COEh  hC  1,2)  +XKMP1,  J,  1)*XI<  I,  J,  2) 

1 +XICMP1, J, 2)*XI< I, J, 1) 

CONTINUE 

COEhK  I,  1 ) = COEh K I , 1)/Ch 
COEhK  I,2)  = C0EhKI,2)/Ch 
CONTINUE 

TYPE  619  5,  C(  COEhK  I,  J),  J»  1,2),  1 = 1,  NUMANT) 


SI  GO  14. BOH 


f! 6-0 Cl  - 7 V 1 3 : 4K  : 14  PACK  ] 6 


a>19S  >OKMAT(  ’ COLB  B * , /*  ?<  1-20.  lb)) 
C 

HETUHM 

HMD 
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■ 'i 


[i 

E 1 

, 


i 
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■ i 


C 

C 

c 


c 


D 

D1000 

C 


20 

C 


40 

C 


2010 

2020 


2030 


C 


C 

D 

UB074 

C 


207  0 
70 
50 
C 


PHIS  PACE  IS  BEST  QUALITY  FIUCIICABU 

SUBROUTINE  SIGOPT  fBQI  OOPI  PUR1USHKD  TO  DDC  ^ 

REPORT  THE  SIMULATION  RESULTS 

IMPLICIT  DOUPLE  PRECISION  CA-H,0-Z) 

INTEGER  CASNUM,  KEI  NI  T,  UPDATE 

COMM  ON/ CON  SNT/ PI  , TWOPI  , ONEOPI 

COMMON /ANSWER/  COMPUTC  10,  1 0)  , I COM  PI  , I C0MP2 

COMMON /SAVE/  DIRVLCC  16,  16,2),XMC  16,  16,  2),  SI  UVARC  16),  DELI- AC 
1 , SI  GLOC(  16) 

C0MM0N/H2H/CASNUM,  I RND,  UKND,  HEIN  I T,  UPDATE,  I TTY,  NUMANT,  I LOOP 
1,NUMSI  G, MSET 

DIMENSION  IORDC  16,  17),  I RPC  16),ILRC  16) 

TYPE  1O0O,  CC  COMPUTC  I, J), 1-1,1  C0MP2)  , J=  1,1  COMF1 ) 
hORMATC  5C  1PD14.  7>) 

DO  20  1=  1,NUMSI  G 

SI  GLOCC  I )*SI  GLOCC  I )*ONEOPI 

CONTINUE 

CALL  SORTDCSIGLOC,NUMSIG,  I OK  DC  1,  1 7)  , I RD,  IL13) 

DO  40  I = 1 , I C OM P 1 

CALL  SORTDC  COMPUTC  1,1  ) ,NUMSI  G,  I OHDC  1,  I ) , I RD,  I LB) 

CONTINUE 

DO  50  1=1 » NUMSI G 
WRI  TEC  2,  2010)1 

10RMATC ///,  * SIGNAL  SOURCE  ’,12,  ' PARAMETERS  :',/) 

WRI  TEC  2,  2020)  SI  GLOCC  I ORIX  I , 17)), SI UVARC  IOKLX  1,17)) 
fORMATC  ’ ANGULAR  LOCATI  ON  =',H('.7i'  * PI  RADIANS', 

1/,'  V/ARI  AN CE  C POWER)  = ' , 1 PD1  0.  3, // ) 

WRI  TEC  2,2030) 

EOhMATC  ' COMPUTED  VALUES  1-OR  ANGULAR  LOCATION  :',// 

1,  ' SAMPLES  THETA/ PI  ERROR/PI',/ 

2,  w RADIANS  RADIANS',/) 


NUMTMS-REINI  T/UPDATE 
N SAM-0 

DO  70  J-  1,NUMTMS 
N SAM- N SAM ♦UP DATE 
ACC-  COMPUTC  I , I ORDC  J,  I ) ) - SI 
TYPE  2074,  IORDC  J,  I ),  I ORDC  I 
I-ORMATC  ' I ORDC  J,  I)-',I3,  ', 
ACC-  COMPUTC  I , I ORDC  I , J)  ) - Si 
ACC- COMPUTC  IORDC  I,J),  J)-SI 
WRI  TEC  2,  2070) NSAM,  COMPUTC  I 
1-ORMATC  2X,  I 5,  10X,H0.7,  6X, 
CONTINUE 
CONTINUE 


GLOCC  IORDC  1,17)) 
, J) 

I ORDC  I , J)  - ',  I 3) 
GLOCC  IORDC  1,17)) 
GLOCC  I ORDC  I,  1?)> 
ORDC  I , J),  J),  ACC 
Me. 7) 


RETURN 

END 

SUBROUTINE  SOITDCX,M,  I ORD,  I H13,  I LB) 
DOUBLE  PRECISION  X 


L 


J. 


' 1 
1 


SI  U01  b.  FOR 
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THIS  PA8S  IS  BEST  QUALITY  mCUCAWil 
FROM  00*"*  FUWUSHEU  TO  DOG  — 

DIMENSI  ON  X ( M > * I OFiDCM)*  I hP(«>,ILWM) 

D TYPE  100*X 

DleO  FORMAT!  ’ X«’,3H0.7) 

C 

IFCM.GT. 1)  GOTO  1 
I OF.  DC  1 > = 1 
RETURN 

1 CONTINUE 
ILK  i>-0 
IhW  1 ) ■>  0 

DO  1A  I = 2,M 
I L PC  I ) *=  0 
IhBC  I )=0 

2 J-  1 

A IF(XCI).GT.XCJ>)  GOTO  10 

6 IFCILBC  JJ.EO.O)  GOTO  H 

J«=  I LBC  J) 

GOTO  A 

8 IKWD-J 

ILKJ)  = I 
GOTO  1 A 

10  I HI  R13C  i)).Lh*C)  GOTO  12 

J=  I RHC  J) 

GOTO  A 

12  I RBC  I > *=i  RH  j) 

I RDC  J)  a I 
1A  CONTINUE 

L=  1 

16  J=  1 

GOTO  20 
18  J-ILPCJ) 

20  I F<  ILBC  J)  .GT.0)  GOTO  1R 

22  IORDCL)»J 

L“L+  1 

2A  IKIRWJ))  28*30*26 

26  JelRBCJ) 

GOTO  20 

33  J--IRBCJ) 

GOTO  22 

30  CONTI  NUE 

D TYPE  31, X 

E31  FORMATC  ' XS«'*3F10.7) 

RETURN  “ 

END 


SI  GO  1 6.  K)R 
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SUBROUTINE  SIG016 
C 

C ABOVE  THRESHOLD  EIGEN  SiL,  ECTI  ON 

C 

IMPLICIT  DOUELE  PRECI SI ONC A-H, O-Z ) 

INTEGER  P.EINIT,  UPDATE,  CASNUM 
COMM  ON /COMCOE/  COEFFC  16,  S),THhESH 
COMMON /JACO/XM U(  32,  32),  El  GEN(  32) 

COMMON/  CONSNT  / PI  , TVOPI  , ON EOPI 

COMM  ON/ Ol  6/  AINC  16,  2),S1GC  1 6,  2)  , SNC  16,2) 

1,  VC  16,  16,2) 

C0MM0N/H2H/CASNUM,  I RN  D,  JRND,  REI  NI  T,  UPDATE,  I TTY,  NbMANT,  I LOOP 
1 , NUMSI  G,MSET, METHOD 
C 
C 

DIMENSI  ON  VC  1 6,  1 6,  2) 

C 

C 

C SELECT  EIGEN  VALUES  ABOVE  THRESHOLD 

C LET  THE  NUMBER  OF  SETS  BE  MSET 

C 

3)99  I FC  THRESH.  NE.  0.  0D0)  GOTO  4000 
TYPE  4010 

4C10  EORMATC/,  ’ THRESHOLD  F OR  SELECTING  EIGENVECTORS  ? ',  S) 

ACCEPT  4015,  THRESH 
4015  F ORMATC  F 1 0.  0) 

GOTO  4020 
4000  CONTINUE 

type  61210,  thresh 

61210  1 ORMATC  /,  ' THRESHOLD  FOR  SELECTING  EIGENVECTORS  :',lPDg.l> 

4020  CONTINUE  “ 

C 

C PASS  1 DATA 

C 

NMANT2=NUMANT 
NMANT1*NMANT2- 1 
I FC METHOD.  EQ.  1 ) NMANT1  = NUMANT 
C TYPE  1 1,NMANT1 

Cll  EORMATC  * NMANTIs  ',  I 4) 

C 

MSETI=0 

DO  6100  M= 1, 2+NMANT1, 2 
I F(  El  GENCM)  . GE.  THRESH)  MSET1=MSET1+ 1 
6100  CONTINUE 

I FCMSETl.  GT.  0)  GOTO  6120 
C 

TYPE  6110,  THRESH,  < El  GENC  I ) , I = 1 , 2*NMANT  1 ) 

6110  EORMATC'  NO  EIGEN  VALUES  ABOVE  THRESH  OF  S1PD20.13, 

1,/,'  EI  GEN  VALUES  ARE' , /,  4C  D20.  1 3)  ) 

STOP*'  THRESH  ERROR* 

C * 

6120  CONTINUE 

C TYPE  6121  1,MSET1 

C61211  FORMAK///,'  ESTIMATED  NUMBER  OF  SOURCES  »',I3) 

C TYPE  6121, THRESH, MSET1,NMANT1 * 


SI  GO!  6.  FOR 
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THIS  Pi 05  IS  BEST  QUALITY  rftAJIl JAJ44 

m\*  dopy  pukhished  ro  duo  

C6121  FORMAT!  * SIG016-  THRESH®  ' , L'l  S. 8 , ' MSL1-M6.  '#  NMAMT  1 ■ 1 1>) 
1)0  6\ SO  1*0, NM AM T I 
K»  0 

RO  61  SO  J«  1, 8*  M S ET 1 , 8 
K-K*  1 

V!  I,K,  1>*XMUC  I,  J) 

VC  I , K,  2)  =XMUC  I 4-NMANT1,  J) 

61  SO  CONTINUE 
617  0 CONTINUE 

C TYPE  6670#  ! ! ! VC  I / J,  K)  , K»  1 * 8)  , Ju  1 * MSET 1 ) , I ® 1 , NMANT 1 > 

C6670  FORMAT!  * V«  ’,  /,  C IX#  ADI  b.R  ) ) 

I FC  METHOD.  EQ.  1 ) GOTO  9000 
C 

C READ  IN  PASS  8 DATA 

C 

C TYPE  8 

FOFiMATC  ' NOW  HEAD  Pi'  DATA’) 

CALL  ASSIGN!  1,  ’ SI  GP2.  TMP  ’ , 9 > 

REA IX  1,  lOJNMANTC 
10  FORMAT!  I 3> 

REAL'!  1,  1 ) ! ! XMU!  I , J> , Ju  1 , 8*NMAN TO)  , I » 1 , 2*NM ANT25 
1 FORMAT! IX, ADI S.8) 

REA IX  1,  1)!  El  GEN!  I ),  I*  1,2*NMANTP> 

CALL  CLOSE!  1) 

C TYPE  3 

C3  FORMAT!  * F I N I SHED  READ  * ) 

C TYPE  A,  < KI  GEN!  I ) , I * 1 , 2*NMAN  Tfi) 

CA  FORMAT!  ’ El  GLN2  ’ , /,  AF  1 0.  7) 

C 

C SELECT  PASS  TWO  DATA  GREATER  THAN  THRESH 

C 

MSET2-0 

DO  7100  M*=  I,  P*NMANT2,  2 
I F!  El  GF'N!  M)  • GE.  THRESH)  MS  ETC*  MS  ETC*  1 
7100  CONTINUE 

I F(  MSET8*  GT.  0)  GOTO  7180 

C 

TYPE  7110,  THRESH,!  El  GfcNC  I ),  I “1,2* NNANT8) 

7110  FORMAT!’  NO  EIGEN  VALUES  AROVh  THRESH  OF  ’,1PDP0.13, 

1,/,’  El  GEN  VALUES  ARE’, /,  A!  D80.  I 3)  ) 

STOP'*  THRESH  ERROR’ 

C 

7120  CONTINUE 

C TYPE  7121  1 , MSF.T2 

C71211  FORMAT!///, ’ ESTIMATED  NUM HER  OF  SOURCES  :’»I3> 

C TYPE  7121,  THRESH,  MSET2,  NMANT2 

C7121  FORMAT!  ’ SIG016-  THRESH®  ’*  D1  5*H,  ’ NSET8*  *,  I 6,  ’,  NMANT2-  ’,  I b> 
DO  /ISO  I-  WNMANT2 
K-0 

DO  7 ISO  J»  1,  £*MSET2#  8 
K-K4  1 

Vl(  1 , K,  1 ) “XMU!  1 » J) 

VIC  l,  K,  2)“XMU!  I 4NMANT2*  J) 

71  SO  CONTINUE 

7170  CONTINUE 


SI  GO  1 6*  h OR 


I 
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OTIS  PACE  IS  BEST  QUALITY  PRACIICABLJ 

C ' K40JI OOPY  /URAI3HKD  TO  HOC  — ^ 

-C 

C 

I VCOL=  0 

DO  2000  J2=2,  2+MSET2+  1 
1 VS=  1 

I FCMODC  J2,  2)  .ME.  0>  GOTO  500  ’RETAIN  COLUMN  TWICE 

I WCOL=I WCOL+ 1 
I WS=0 

500  CONTINUE 

CALL  ORGC  W,  NMANT 1 » M SET 1 » I V/COL,  I WS) 

2000  CONTINUE 

C TYPE  1020, CCC  VCI,  J,K),K=1,2),  J= 1 , MS ET 1 ) , I = 1 , NMANT 1 ) 

Cl 020  FORMAT!  ’ RESULANT  V= ’, /, C IX, ADI  5.8 ) ) 

NUMANT=NMANT1 
9000  MSET=MSET1 

TYPE  8 121  1,MSET 

81211  FORMAT!///,  * ESTIMATED  NUMBER  OE  SOURCES  :*,I3) 

CALL  SIGABT* 

RETURN 

END 

SUBROUTINE  ORGC  V,  NM AN T,  MSETV,  I WCOL,  I VS) 

C 

c 

c 

IMPLICIT  DOUBLE  PRECISION  CA-H,0-Z) 

COMMON/SAVE/DI  KVECC  16,  16,  2>,XM(  16,  1 6,  2)  , SI  GVARC  16),  DELI- AC 
1,  SI  GLOCC  1 6),  VLIMI  T 

COMMON/Ol  6/AC  16,  2),BC  16,  2),SN<  16,  2),  VC  16,  16,2) 

C 

DIMENSION  WPC  16,  2),  VC  16,  16,  2) 

C 

C TYPE  9 00,  C C VC  I , 1 , K ) , K=  1,2), 1=1,  NMANT) 

©00  FORMAT!  ' ORG-VC  1 TO  NMANT,  1,1  TO  2 ’, /AC  IX,  1 PD1  5.8  ) > 

C TYPE  910,  C(W(I,I  VICOL,  K>  , K=  1,2),  1=  1+1  WS,  NMANT+ 1 V/S) 

©10  FORMAT!  * ORG-W',/,  A!  IX,  1PD15.8)) 

C XR1  = 0.0&0 

C X R2=  0*0  D0 

C DO  920  I * 1,  NMANT 

C DO  920  1, MSETV 

C XR1=XR1  + V!  I,  J,  1) 

C XR2«=XR2+ VC  I , J,  2) 

©20  CONTINUE 

C TYPE  925,XR1,XR2 

©25  FORMAT!  ’ ORG-SUM  OF  V=  ' , 2!  IX,  1PD1  5.8 > ) 

DO  1 000 " J= 1, MSETV 
A!  J,  1 ) = 0 » 0 DO 
A!  J,  2)  = 0.  0D0 
DO  100  I - 1, NMANT 

CALL  CMPLXM!  Vi!  I +1  WS,  I WlCOL,  1),  W!  I+I  V/S,  I WCOL,  2),  V!  I , J,  1) 

1,-VC  I,  J,  2),XR1,XR2> 

AC  J,  1 ) s AC  J,  1>+XR1 
AC  J,  2>*AC  J,  2)  +XR2 
100  CONTINUE 

1000  CONTINUE 


SI  GOl  6.  FOR 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABI4 
FROM  OGPY  JUfUUSHED  TOJ3DO  — 

V 

C SET  B = SUM  ACJJVCJ) 

C 

DO  50  I = 1#NMANT 
EC  I,  1>=0.0D0 
B<  I<  2>  = 0.0D0 
DO  50  J=  1#  MSETV 

CALL  CMPLXMC  AC  J#  1 ) , AC  J#  2)  # VC  I * J#  1 ) # VC  I , J#  2>#XR1#XR2) 

BC  I*  1>  = BC  I#  1)+XR1 
BCI#2>  = BCI#2)+XR2 
50  CONTINUE 

C 

C WP=W-SUMOFC  BCJJVCJ)] 

c 

DO  70  1 = 1#  NMANT 

wpc  i # n = wc  i + i ws#  i wcol#  n-Bc  i » n 
70  Wpc  I # 2)  =wc  I + 1 WS#  I WOOL#  2)  - BC  I # 2) 

C TYPE  30  5#  C C WPC  I #J)#J=  1#  2>#  1=  1#NMANT) 

C305  FORMATC  ’ VP=  • , /#  C ADI  2*  5)  ) 

C " 

C INNER  PRODUCT 

C 

C DO  3000  I = 1#MSETV 

C AC  I#  1 ) = 0»  0D0 

C ACI#2)  = 0. 0D0 

C DO  3000  J=l# NMANT 

C CALL  CMPLXMC  WPC  J#  1>#-WPC  J#2)#VC  J#I#  1),  VC  J#I#2)#XR1#XR2> 

C AC  I#  1)=ACI#  1)  +X R 1 

C AC  I#2>  = AC  I#2>  +XR2 

C3000  CONTINUE 

C TYPE  3010#CCACI#K)#  K=  1 # 2)  # I = 1 #M  SETV) 

C3010  KORMATC*  INNER  PRODUCT  OF’  VP*  V=  ’ # /,  4C  1 PD1  5.8  > > 

C 

VX  1=0.  0 

DO  400  1 = 1, NMANT 

VX1  = VX1  + VPCI#  1)*WPCI#  1)  + WPC  I#2)*VPC  1,2) 

400  CONTINUE 

1=  1 

IFCVXI.LT.WLIMIT)  1 = 0 
D TYPE  405#  WX  1 # I 

D405  FORMATC*  INNER  PRODUCT=  ’ # G 1 0.  2#  ’ # I=’#I2) 

I FC  I .NE.'  1>  GOTO  2000 
C 

C MOVE  WP  INTO  V 

C 

MSETV-MSETV+1 

ODSVX  1=1*  0D0/DSQRTC  WX  1 ) 

DO  8 0 J=  1# NMANT 
VC  J#  MSETV,  1)=WPC  J#  1 ) * ODSWX  1 
VC  J#  MSETV#  2)  = WPC  J#  2)  * ODSVX  1 
80  CONTINUE 

C 

2000  CONTINUE 

RETURN 
END 


f 
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SUFHOUTI  M E SI  GO  17 
METHOD  TWO 

1MFHCIT  DOUBLE  PRECISION  CA-H,0-7.) 

INTEGR.R  CASNUM,  RU  .\J  1 T,  UPDATE 

COMM  ON/COMCOE/ COER  K 1 6,  O) , THRESH  * AM  TNI  S 

COMMON /HKH/ CASNUM,  I RN  D»  JRND,  KR.l  N1  T»  UPDATE,  I TTY  » N l WANT 

1 , I L OOP,  MUM  SI  G,  MS  ET,  M ETH  0 D 

COMMON /JACO/XMU<  32,  32)  , El  GRSJC  32) 

DIMENSION  XMU2C  32,  32),  El  GEN2C  32)  ,AC  16,C),A2C  16,2) 

DIMENSI  ON  V(  1 6,  1 6,  2)  , VC  1 6,  2)  , ASC  1 6)  , AS2C  1 6) 


READ  PASS  TWO  DATA Vi 

CALL  ASSIGN!  1,  ’ SI  GP8»  TMl'  * , 9 ) 

READ!  1,  1 0)  NM  AN  Til 
RORMATC  I 3) 

KR  ADC  1 , 1 > C C XM02C  I , J)  , J<=  1 , 2*NMAN T2)  , 1 - 1 , B*NMAN  T2) 

RORMATC  IX,  ADI  S-H  ) 

lit  AL-C  1 , 1 ) C El  GEN2C  I ) , I - 1 , 2*tNMANTfi) 

CALL  CL  OS  EC  1) 

NM  AN  T 1 « N M AN  T 2-  1 
DO  61  SO  I*>l,NMANTl 
K-  0 

DO  6 ISO  1, 2*NMANT1, 2 

K-K>  1 

VC  I , K,  1)-XMU(  I,J) 

V(  I , K • £>-XMUC  I +NMANT1,  J) 

CONTINUE 

DO  7 ISO  J«  1 , NMAN  T2 
WC  J,  1)-XMU2C  J,  1) 

WC  J,  8)  “XMU2C  J*NMANT2,  1) 

CONTINUE 

INNER  FHODUCT 

DO  100  J- 1 , NMANT I 
AC  J,  l)-0.  ODD 
AC  J,  8)«0.0D0 
ABC  J,  1 > “ 0«  ODO 
A8C  J,  2)  • 0»  ODO 
DO  100  I NMANT  1 

CALL  CMI'LXMC  VC  1 , 1 ) , - Vl  I , 2)  , VC  I , J,  1 ) , VC  I , J,  2)  , XK 1 , XR2) 

AC  J,  1 ) - AC  0,  1 )+XHI 
AC  J,  8)  - AC  J,  8)+XR2 

CALL  CMPLXMC  ViC  I ♦ 1,  1 ) , - WC  I ♦ 1 , 2)  , VC  I , J,  1 ) , VC  I , J,  2)  , XR1  ,XR2) 
A2C  J,  1 ) - ABC  J,  1 ) +XR1 
ABC  J,  2)- ABC  J,  8) ♦X H2 
CONTINUE 

TYPE  101»C  C AC  I ,K),K*  1, 8),  I - 1 ,NMANT1 ) 
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D'101  FOFiMATC  • INNER  PRODUCT  AT  VECTOR  1 /,  C IX#  4D1  5.8  ) ) 

D TYPE  1B2,CC  A2C  I , K ) , K=  1 , 2)  , I = 1 , NMAN T 1 > 

D102  FORMATC  ’ INNER  PRODUCT  AT  VECTOR  2 ’ , /,  C IX#  4D1  5*8  ) ) 

DO  8 100'l  = 1,  NMANT  1 

ASC  I ) = DSGRTC  AC  I * 1)  i AC  I , 1 )+AC  I * 2)  * AC  1,2)) 

AS2C  I ) = DSGRTC  A2C  I , 1 )*A2C  I , 1 )+A2C  I , 2)*A2C  1,2)) 

8160  CONTINUE 

TYPE  8 IBS,  C ASC  I ),  1 = 1,  NMANT  1) 

8165  FORMATC/,  * INNER  PRODUCTS  WITH  TEST  VECTOR  1 : /4C  IX,  D1  5*8  ) > 

TYPE  8 106,'C  AS2C  I ),  i = ;,nmantd 

8166  F ORMATC  /,  * INNER  PRODUCTS  WITH  TEST  VECTOR  2 : ’, /4C  IX,  D1  5. 8)  ) 
TYPE  3610' 

3010  F ORNATC  //,  * NUMBER  OF  EIGENVECTORS  TO  BE  SELECTED  EASED  ON*, 

1,  * INNER  PRODUCTS  ? ’) 

ACCEPT  301 5, MSET 
3015  F ORMATC  I 2) 

C 

6120  CONTINUE 

DO  6170  1=1, NMANT 1-1 

DO  6170  J= 1, MSET 
VCI,J,  2)=-VCI,  J,  2) 

617  0 CONTINUE 

SUM SGW=  0.  0D0 
I=NMANT1 

DO  6180  J=1,MSET 

SOMSQW=SUMSQW+ VC  I , J,  1 ) * VC  I , J,  1)+VC  I , J,  2)  * VC  I , J,  2) 

618  0 CONTINUE 

D TYPE  618  1,  SUMSQW 

D618  1 F ORMATC  ’ SUMSCW=  ’ , 1 PD20.  1 2) 

I FC  DAPSC  1 . 0D0-SUMSGW)  .LT.  1 . 6D-5)  TYPE  618  5 
6185  F ORMATC  /,  * WARNING  ABSOLUTE  VALUE  OF  I-SUMSQV  IS  LESS  THAN  1 •CD- 
IS’) 

C 

C COMPUTE  COEFF 

C 

COEFFC  NMANT  1,  1)  = 1.0D0 
CF=-C 1 • 0D0-SUMSQW) 

COEFFC  NMANT  1, 2)  = 0.  0D0 

MP1=NMANT1 

DO  6199  I = 1, NMANT l-  1 

COEFFC  I,  1)  = 0.  0D0 

COEFFC I, 2)=0. 0D0 

DO  6190  J=l, MSET 

COEFFC  I,  1)  = C0EFFCI,  1)+VCMP1,J,  1 ) * VC  I , J,  1) 

1 ‘ -VCMP1,J,2)*VCI,J,2) 

COEF  FC  I , 2)  = COEFFC  I,2)+VCMP1,J,  1)*VC  I,  J,  2) 

1 +VCMP1,  J,  2)*VC  I,  J,  1) 

619  0 CONTINUE 

COEFFC I, 1 ) =COEFFC I , 1)/CF 
C0EFFCI,2)  = C0EFFCI,2)/CF 
6199  CONTINUE 

C 

D TYPE  6195,CC  COEFFC  I,  J),  J=  1 , 2)  , I = 1 , NMANT  1 ) 

D6195  F ORMATC  ’ COEFF’,/,  2CF20.  15)) 

C * 


/ 


SIG017. FOR 


26- OCT- 7 7 14:10:2b  PAGE 


25 


NUMANT=NMAMT1 

RETURM 

END 
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SUBROUTINE  SIG018 

ABOVE  THRESHOLD  El  GEN  SELECTION 


THIS  PACE  IS  BEST  QUALITY  PEACUCA&H 
PROA  COPY  PUKAISHSU  TV DOC  - 


4010 

4015 

400  0 

61210 

4020 


6100 


6110 


C 

6120 

6121 1 
D 

D6121 


6150 

C 


IMPLICIT  DOUBLE  PRLCI  SI  ON(  A-H  > O-Z) 

INTEGER  REINIT,  UPDATE,  CASNUM 
COMMON/COMCOL/  COEF  PC  1 6,  2)  , THRESH#  ANTNI  S 
COMMON / JACO/XMU!  32,32),  LI  GENC  32) 

COMMON/  CONSNT  / PI  , TVOPI  , ONEOPI 
DIMENSION  XI  <16,  16,  2) 

COMM  ON /H  2.4/ CASNUM,  I END,  JRN  D,  REI  NI  T,  UPDATE,  I TTY,  NUM  ANT,  I LOOP 
1, NUMSI  G,MSET,METHOD 

SELECT  EIGEN  VALUES  ABOVE  THRESHOLD 
LET  THE  NUMBER  OP  SETS  BE  MSET 

I K THRESH. NE.  0.  0D0)  GOTO  4000 
TYPE  4010 

FORMAT!/,  * THRESHOLD  FOR  SELECTING  EIGENVECTORS  ? *,  1) 

ACCEPT  401  5,  THRESH 
FORMAT!  F10.  0) 

GOTO  4020 
CONTINUE 

TYPE  61210,  THRESH 

FORMAT!/,  * THRESHOLD  FOR  SELECTING  EIGENVECTORS  1PD8.1) 
CONTINUE  " 

MSET=0 

DO  6100  M=  1, 2+NUMANT,  2 

I F<  El  GEN! M)  . GE.  THRESH)  MSET=MSET+  1 

CONTINUE 

I F<MSET.GT.0)  GOTO  6120 

TYPE  6110,  THRESH,!  El  GEN!  I ),  1 = 1, 2*NUMANT) 

FORMAT!  * NO  EIGEN  VALUES  ABOVE  THRESH  OF  *,1PD20.13, 

1, /,  * EfGEN  VALUES  ARE*, /,  4!  D20.  1 3)  ) 

STOP"’  THRESH  ERROR* 

CONTINUE 
TYPE  61211, MSET 

FORMAT!///,  * ESTIMATED  NUMBER  OF  SOURCES  :*,I3) 

type  6121,  tf?resh,mset,nimant 

FORMAT!  ' SIGABT-  THRESH=  *,D15«8,  * MSET=*,I6,  ’,  NUMANT-  *,  I 5) 
DO  6150 '1  = 1, NUMANT 
K=  0 

DO  6150  J=1,2*MSET»2 
K=K+  1 

XI < I , K,  1)=XMU!I, J) 

XI!  I,K,2)=XMU<  I+NUMANT,  J) 

CONTINUE 

DO  6170  1 = 1 , NUMANT-  1 
DO  6170  J= 1, MSET 
XI!I, J, 2)=-XI!I,  J, 2) 

CONTINUE 


6170 

C 
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C 

0 


618  0 
D 

D618  l 

618  5 

C 

C 

C 


619  0 


6199 

C 

D 

D619S 

C 


COMPUTE  LAMBDA* 

SUMSQW=  0. 0D0 
I-NUMANT 

DO  6180  1 * MSET 


27 
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SUM  OF  THE  SQUARES  OF  Vi 


SUMSQV=SUMSQV+XI( I, J,  1 > *XI  C I * J*  1 > +XI  < I * J*  2)  *X I C I * J*  2) 

CONTINUE 

TYPE  618  1*  SUMSQW 

FORMAT!  ' SUMSQW*  '*  1PD20.  12) 

IKDABSf  1.0D0-SU1^SQW).LT.  1.0D-5)  TYPE  6135 

FORMAT!/,’  WARNING  ABSOLUTE  VALUE  Of  1-SUMSQV  IS  LESS  THAN  l.flD- 
15*) 


COMPUTE  COEFF 


COEFF!  NUM  ANT*  1)«=1.0D0 
CF  = - ! 1.0D0-SUMSQW) 

COEFF!  NUMANT*  2)  = 0*  0D0 

MP1«=NUMANT 

DO  6199  I = 1*NUMANT-1 

COEF  F!  I,  1)  = 0.0D0 

COEFF!  I * 2)  = 0.  GD0 

DO  6190  J=  1 * M SET 

COEF  FC  I * l)  = COEFF!  I , 1 ) +XI ! MP1  * J*  1>*XI!I*J*  1) 

1 -XICMP1*  J*2)+XI(I*J*2) 

COEFF!  I*2>  = COEFF!  I * 2)  +XI  ! MPl » J*  l)*XI! I*  J*  2) 

1 +XI!MP1*  J,2>*XI!  I*  J*  1) 

CONTINUE 

COEFF!  I,  1 ) = COEFF!  I * 1)/CF 
COEF  F(  I , 2>  = C0EFF!  1*2)  /CF 
CONTINUE 

TYPE  619  5*C(C0EFF(I*J),J=l*2)*I  = l,NUMANT) 
FORMAT!  ’ COEFF'*/*  2(F20.15)) 

RETURN 

END 


